diff --git a/docs/cpp_api/continuum_mechanics_overview.md b/docs/cpp_api/continuum_mechanics_overview.md index dded701ed..5db1aea97 100644 --- a/docs/cpp_api/continuum_mechanics_overview.md +++ b/docs/cpp_api/continuum_mechanics_overview.md @@ -80,8 +80,7 @@ User Material subroutines (UMAT) for finite element analysis, organized by strai - `damage_weibull` - Weibull statistical damage model **Shape Memory Alloys (SMA):** -- `SMA_mono` - Monolithic SMA model -- `SMA_mono_cubic` - Cubic crystal SMA model +- `SMA_mono` - Micromechanical monocrystal SMA model (supports SMAMO isotropic, SMAMC cubic symmetry) - `unified_T` - Unified thermomechanical SMA model **Combined:** diff --git a/environment_doc.yml b/environment_doc.yml index 2c01837f7..93895b575 100644 --- a/environment_doc.yml +++ b/environment_doc.yml @@ -38,3 +38,4 @@ dependencies: - pandas - breathe - doxygen + - scikit-learn \ No newline at end of file diff --git a/examples/analysis/hyperelastic_parameter_identification.py b/examples/analysis/hyperelastic_parameter_identification.py new file mode 100644 index 000000000..82ee0b723 --- /dev/null +++ b/examples/analysis/hyperelastic_parameter_identification.py @@ -0,0 +1,828 @@ +""" +Parameter Identification for Mooney-Rivlin Hyperelastic Material +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This example demonstrates inverse parameter identification for hyperelastic +materials using experimental stress-strain data. We identify the Mooney-Rivlin +model parameters from Treloar's classical rubber elasticity experiments. + +**Optimization**: scipy.optimize.differential_evolution (global optimizer) +**Cost function**: Mean Squared Error on stress (sklearn.metrics.mean_squared_error) +**Forward model**: simcoon hyperelastic stress functions +""" + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from scipy.optimize import differential_evolution, fsolve +from sklearn.metrics import mean_squared_error +from simcoon import simmit as sim +import os + +################################################################################### +# Methodology Overview +# -------------------- +# +# **Inverse Problem Formulation** +# +# Parameter identification is an inverse problem: given experimental observations +# (stress-strain data), we seek the material parameters that minimize the +# discrepancy between model predictions and experiments. +# +# The optimization problem is: +# +# .. math:: +# \boldsymbol{\theta}^* = \arg\min_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) +# +# where :math:`\boldsymbol{\theta}` are the material parameters and +# :math:`\mathcal{L}` is the cost function. +# +# **Cost Function: Stress-Based MSE** +# +# We use the Mean Squared Error (MSE) computed on the first Piola-Kirchhoff +# stress :math:`P_{11}`: +# +# .. math:: +# \text{MSE} = \frac{1}{N} \sum_{i=1}^{N} \left( P_{11}^{\text{model}}(\lambda_i; \boldsymbol{\theta}) - P_{11}^{\text{exp}}(\lambda_i) \right)^2 +# +# This stress-based metric directly measures the model's ability to predict +# mechanical response, which is the primary quantity of interest in constitutive +# modeling. +# +# **Why Differential Evolution?** +# +# We use ``scipy.optimize.differential_evolution`` because: +# +# 1. **Global optimization**: Unlike gradient-based methods, it explores the +# entire parameter space and avoids local minima +# 2. **Derivative-free**: No gradients required, robust for non-smooth objectives +# 3. **Bounded search**: Natural handling of physical parameter constraints +# 4. **Robustness**: Less sensitive to initial guess compared to local methods +# +# For hyperelastic models, the cost function landscape can be non-convex, +# making global optimization essential for reliable parameter identification. + +################################################################################### +# Mooney-Rivlin Constitutive Model +# -------------------------------- +# +# The Mooney-Rivlin model is a phenomenological hyperelastic model widely used +# for rubber-like materials. The strain energy function is: +# +# .. math:: +# W = C_{10}(\bar{I}_1 - 3) + C_{01}(\bar{I}_2 - 3) + U(J) +# +# where: +# +# - :math:`\bar{I}_1, \bar{I}_2` are the first and second isochoric invariants +# of the left Cauchy-Green tensor :math:`\mathbf{b} = \mathbf{F}\mathbf{F}^T` +# - :math:`C_{10}, C_{01}` are the material parameters to identify +# - :math:`U(J) = \kappa(J \ln J - J + 1)` is the volumetric contribution +# - :math:`\kappa` is the bulk modulus (fixed, assuming near-incompressibility) +# +# The Cauchy stress is computed as: +# +# .. math:: +# \boldsymbol{\sigma} = \boldsymbol{\sigma}_{\text{iso}} + \boldsymbol{\sigma}_{\text{vol}} +# +# where simcoon provides ``sigma_iso_hyper_invariants`` and ``sigma_vol_hyper`` +# to compute these contributions from the strain energy derivatives: +# +# .. math:: +# \frac{\partial W}{\partial \bar{I}_1} = C_{10}, \quad +# \frac{\partial W}{\partial \bar{I}_2} = C_{01}, \quad +# \frac{\partial U}{\partial J} = \kappa \ln J + +################################################################################### +# Experimental Data: Treloar's Rubber Experiments +# ----------------------------------------------- +# +# L.R.G. Treloar (1944) performed classical experiments on vulcanized rubber +# under three loading conditions: +# +# 1. **Uniaxial Tension (UT)**: :math:`\lambda_1 = \lambda, \lambda_2 = \lambda_3 = \lambda_t` +# 2. **Pure Shear (PS)**: :math:`\lambda_1 = \lambda, \lambda_2 = \lambda_t, \lambda_3 = 1` +# 3. **Equibiaxial Tension (ET)**: :math:`\lambda_1 = \lambda_2 = \lambda, \lambda_3 = \lambda_t` +# +# where :math:`\lambda_t` is the transverse stretch determined by equilibrium +# (zero transverse stress). The data provides stretch :math:`\lambda` and the +# corresponding first Piola-Kirchhoff stress :math:`P_{11}` [MPa]. + + +def load_treloar_data(filepath: str) -> pd.DataFrame: + """ + Load Treloar experimental data. + + The file contains columns for three loading cases: + - lambda_1, P1_MPa: Uniaxial tension + - lambda_2, P2_MPa: Pure shear + - lambda_3, P3_MPa: Equibiaxial tension + + Parameters + ---------- + filepath : str + Path to Treloar.txt data file + + Returns + ------- + pd.DataFrame + Experimental data with columns for each loading case + """ + df = pd.read_csv( + filepath, + sep=r"\s+", + engine="python", + names=["lambda_1", "P1_MPa", "lambda_2", "P2_MPa", "lambda_3", "P3_MPa"], + header=0, + ) + return df + + +################################################################################### +# Forward Model: Stress Prediction using simcoon +# ---------------------------------------------- +# +# The forward model computes the first Piola-Kirchhoff stress for given +# material parameters and stretch values. This is the core function that +# will be called repeatedly during optimization. +# +# **Key steps:** +# +# 1. For each stretch :math:`\lambda_1`, solve for the transverse stretch +# :math:`\lambda_t` that satisfies equilibrium (zero transverse stress) +# 2. Construct the deformation gradient :math:`\mathbf{F}` +# 3. Compute the left Cauchy-Green tensor :math:`\mathbf{b} = \mathbf{F}\mathbf{F}^T` +# using ``sim.L_Cauchy_Green(F)`` +# 4. Compute strain energy derivatives for Mooney-Rivlin +# 5. Compute Cauchy stress using ``sim.sigma_iso_hyper_invariants`` and +# ``sim.sigma_vol_hyper`` +# 6. Convert to first Piola-Kirchhoff stress: :math:`P_{11} = \sigma_{11}/\lambda_1` + + +def mooney_rivlin_pk1_stress( + C10: float, + C01: float, + kappa: float, + lambda_array: np.ndarray, + loading_case: str = "UT", +) -> np.ndarray: + """ + Compute first Piola-Kirchhoff stress for Mooney-Rivlin material. + + Uses simcoon's hyperelastic stress functions to compute the Cauchy stress + from the strain energy derivatives, then converts to PK1 stress. + + Parameters + ---------- + C10 : float + First Mooney-Rivlin parameter [MPa] + C01 : float + Second Mooney-Rivlin parameter [MPa] + kappa : float + Bulk modulus for volumetric response [MPa] + lambda_array : np.ndarray + Array of principal stretch values in the loading direction + loading_case : str + Type of loading: "UT" (uniaxial tension), "PS" (pure shear), + or "ET" (equibiaxial tension) + + Returns + ------- + np.ndarray + First Piola-Kirchhoff stress P_11 [MPa] + """ + PK1_stress = [] + lambda_t_guess = 1.0 # Initial guess for transverse stretch + + for lam in lambda_array: + # ----------------------------------------------------------------- + # Step 1: Find transverse stretch satisfying equilibrium + # ----------------------------------------------------------------- + def equilibrium_residual(lambda_t_vec): + """Residual function for finding equilibrium transverse stretch.""" + lt = lambda_t_vec[0] + + # Construct deformation gradient based on loading case + if loading_case == "UT": + # Uniaxial: F = diag(lambda, lambda_t, lambda_t) + F = np.diag([lam, lt, lt]) + J = lam * lt**2 + elif loading_case == "PS": + # Pure shear: F = diag(lambda, lambda_t, 1) + F = np.diag([lam, lt, 1.0]) + J = lam * lt + elif loading_case == "ET": + # Equibiaxial: F = diag(lambda, lambda, lambda_t) + F = np.diag([lam, lam, lt]) + J = lam**2 * lt + else: + raise ValueError(f"Unknown loading case: {loading_case}") + + # Compute left Cauchy-Green tensor using simcoon + b = sim.L_Cauchy_Green(F) + + # Mooney-Rivlin strain energy derivatives + dW_dI1_bar = C10 + dW_dI2_bar = C01 + dU_dJ = kappa * np.log(J) + + # Compute Cauchy stress using simcoon functions + sigma_iso = sim.sigma_iso_hyper_invariants( + float(dW_dI1_bar), float(dW_dI2_bar), b, J, False + ) + sigma_vol = sim.sigma_vol_hyper(dU_dJ, b, J, False) + sigma = sigma_iso + sigma_vol + + # Return stress component that should be zero at equilibrium + if loading_case == "UT": + # sigma_22 = sigma_33 = 0 + return 0.5 * (sigma[1, 1] + sigma[2, 2]) + elif loading_case == "PS": + # sigma_22 = 0 + return sigma[1, 1] + else: # ET + # sigma_33 = 0 + return sigma[2, 2] + + # Solve for equilibrium transverse stretch + lambda_t = fsolve(equilibrium_residual, lambda_t_guess, full_output=False)[0] + lambda_t_guess = lambda_t # Use as starting point for next iteration + + # ----------------------------------------------------------------- + # Step 2: Compute final stress state at equilibrium + # ----------------------------------------------------------------- + if loading_case == "UT": + F = np.diag([lam, lambda_t, lambda_t]) + J = lam * lambda_t**2 + elif loading_case == "PS": + F = np.diag([lam, lambda_t, 1.0]) + J = lam * lambda_t + else: # ET + F = np.diag([lam, lam, lambda_t]) + J = lam**2 * lambda_t + + b = sim.L_Cauchy_Green(F) + + # Strain energy derivatives + dW_dI1_bar = C10 + dW_dI2_bar = C01 + dU_dJ = kappa * np.log(J) + + # Compute Cauchy stress + sigma_iso = sim.sigma_iso_hyper_invariants( + float(dW_dI1_bar), float(dW_dI2_bar), b, J, False + ) + sigma_vol = sim.sigma_vol_hyper(dU_dJ, b, J, False) + sigma = sigma_iso + sigma_vol + + # Convert Cauchy stress to PK1: P = J * sigma * F^{-T} + # For diagonal F: P_11 = sigma_11 / lambda_1 + PK1_stress.append(sigma[0, 0] / lam) + + return np.array(PK1_stress) + + +################################################################################### +# Cost Function: Stress-Based Mean Squared Error +# ---------------------------------------------- +# +# The cost function quantifies the discrepancy between model predictions and +# experimental data. We use MSE computed on the PK1 stress values. +# +# **Normalized MSE (NMSE)**: +# +# When combining multiple loading cases with different stress magnitudes, +# we can normalize by the variance to ensure balanced weighting: +# +# .. math:: +# \text{NMSE} = \frac{\text{MSE}}{\text{Var}(P^{\text{exp}})} +# +# This ensures that loading cases with smaller stress magnitudes are not +# dominated by cases with larger stresses during combined fitting. + + +def cost_function_mse( + params: np.ndarray, + lambda_exp: np.ndarray, + P_exp: np.ndarray, + kappa: float, + loading_case: str, + normalize: bool = False, +) -> float: + """ + Compute MSE cost between model prediction and experimental stress data. + + Parameters + ---------- + params : np.ndarray + Material parameters [C10, C01] + lambda_exp : np.ndarray + Experimental stretch values + P_exp : np.ndarray + Experimental PK1 stress values [MPa] + kappa : float + Fixed bulk modulus [MPa] + loading_case : str + Loading type ("UT", "PS", or "ET") + normalize : bool + If True, divide MSE by variance of experimental data + + Returns + ------- + float + MSE or NMSE value + """ + C10, C01 = params + + try: + # Predict stress using forward model + P_model = mooney_rivlin_pk1_stress(C10, C01, kappa, lambda_exp, loading_case) + + # Compute MSE using sklearn + mse = mean_squared_error(P_exp, P_model) + + if normalize: + variance = np.var(P_exp) + return mse / variance if variance > 1e-12 else mse + return mse + + except Exception: + # Return large cost if computation fails (e.g., numerical issues) + return 1e10 + + +def cost_function_combined( + params: np.ndarray, data_dict: dict, kappa: float, normalize: bool = True +) -> float: + """ + Combined cost function for simultaneous fitting to multiple loading cases. + + Parameters + ---------- + params : np.ndarray + Material parameters [C10, C01] + data_dict : dict + Dictionary mapping loading case names to (lambda, P) tuples + kappa : float + Fixed bulk modulus [MPa] + normalize : bool + If True, use NMSE for balanced weighting across cases + + Returns + ------- + float + Sum of (N)MSE values over all loading cases + """ + total_cost = 0.0 + for case_name, (lambda_exp, P_exp) in data_dict.items(): + cost = cost_function_mse(params, lambda_exp, P_exp, kappa, case_name, normalize) + total_cost += cost + return total_cost + + +################################################################################### +# Parameter Identification using Differential Evolution +# ----------------------------------------------------- +# +# ``scipy.optimize.differential_evolution`` is a stochastic global optimizer +# based on evolutionary algorithms. Key parameters: +# +# - **bounds**: Search space for each parameter +# - **strategy**: Mutation strategy ('best1bin' works well for most problems) +# - **popsize**: Population size (larger = more thorough search) +# - **mutation**: Mutation constant (controls exploration vs exploitation) +# - **recombination**: Crossover probability +# - **polish**: If True, refines result with L-BFGS-B (local optimizer) + + +def identify_parameters( + lambda_exp: np.ndarray, + P_exp: np.ndarray, + kappa: float = 4000.0, + loading_case: str = "UT", + bounds: list = None, + normalize: bool = False, + seed: int = 42, + verbose: bool = True, +) -> dict: + """ + Identify Mooney-Rivlin parameters using differential evolution. + + Parameters + ---------- + lambda_exp : np.ndarray + Experimental stretch values + P_exp : np.ndarray + Experimental PK1 stress values [MPa] + kappa : float + Fixed bulk modulus [MPa] (default: 4000 for near-incompressibility) + loading_case : str + Loading type ("UT", "PS", or "ET") + bounds : list + Parameter bounds [(C10_min, C10_max), (C01_min, C01_max)] + normalize : bool + If True, use normalized MSE + seed : int + Random seed for reproducibility + verbose : bool + Print optimization progress + + Returns + ------- + dict + Optimized parameters and optimization results + """ + if bounds is None: + bounds = [(0.01, 2.0), (-1.0, 1.0)] + + if verbose: + print(f"\n{'=' * 60}") + print(f"Optimizing for {loading_case} loading case") + print(f"{'=' * 60}") + print(f"Parameter bounds: C10 in {bounds[0]}, C01 in {bounds[1]}") + print(f"Bulk modulus (fixed): kappa = {kappa} MPa") + + result = differential_evolution( + cost_function_mse, + bounds=bounds, + args=(lambda_exp, P_exp, kappa, loading_case, normalize), + strategy="best1bin", + maxiter=500, + popsize=15, + tol=1e-8, + mutation=(0.5, 1.0), + recombination=0.7, + seed=seed, + polish=True, # Refine with L-BFGS-B + disp=verbose, + ) + + C10_opt, C01_opt = result.x + + if verbose: + print(f"\nOptimization completed!") + print(f" C10 = {C10_opt:.6f} MPa") + print(f" C01 = {C01_opt:.6f} MPa") + print(f" Final MSE = {result.fun:.6e} MPa^2") + + return { + "C10": C10_opt, + "C01": C01_opt, + "kappa": kappa, + "mse": result.fun, + "success": result.success, + "result": result, + } + + +def identify_parameters_combined( + data_dict: dict, + kappa: float = 4000.0, + bounds: list = None, + normalize: bool = True, + seed: int = 42, + verbose: bool = True, +) -> dict: + """ + Identify parameters by fitting to multiple loading cases simultaneously. + + Using NMSE (normalized MSE) ensures balanced contribution from each + loading case, preventing cases with larger stress magnitudes from + dominating the optimization. + + Parameters + ---------- + data_dict : dict + Dictionary {case_name: (lambda_exp, P_exp)} for each loading case + kappa : float + Fixed bulk modulus [MPa] + bounds : list + Parameter bounds [(C10_min, C10_max), (C01_min, C01_max)] + normalize : bool + If True, use NMSE (recommended for combined fitting) + seed : int + Random seed + verbose : bool + Print progress + + Returns + ------- + dict + Optimized parameters and results + """ + if bounds is None: + bounds = [(0.01, 2.0), (-1.0, 1.0)] + + if verbose: + print(f"\n{'=' * 60}") + print(f"Combined optimization for: {list(data_dict.keys())}") + print(f"{'=' * 60}") + print(f"Parameter bounds: C10 in {bounds[0]}, C01 in {bounds[1]}") + print(f"Bulk modulus (fixed): kappa = {kappa} MPa") + print(f"Using {'NMSE' if normalize else 'MSE'} for cost function") + + result = differential_evolution( + cost_function_combined, + bounds=bounds, + args=(data_dict, kappa, normalize), + strategy="best1bin", + maxiter=500, + popsize=15, + tol=1e-8, + mutation=(0.5, 1.0), + recombination=0.7, + seed=seed, + polish=True, + disp=verbose, + ) + + C10_opt, C01_opt = result.x + + if verbose: + print(f"\nOptimization completed!") + print(f" C10 = {C10_opt:.6f} MPa") + print(f" C01 = {C01_opt:.6f} MPa") + print(f" Final combined cost = {result.fun:.6e}") + + return { + "C10": C10_opt, + "C01": C01_opt, + "kappa": kappa, + "cost": result.fun, + "success": result.success, + "result": result, + } + + +################################################################################### +# Visualization Functions +# ----------------------- + + +def plot_fit_comparison( + lambda_exp: np.ndarray, + P_exp: np.ndarray, + params: dict, + loading_case: str, + title: str = None, + ax=None, +): + """ + Plot experimental data vs model prediction. + + Parameters + ---------- + lambda_exp : np.ndarray + Experimental stretch values + P_exp : np.ndarray + Experimental PK1 stress values [MPa] + params : dict + Dictionary with 'C10', 'C01', 'kappa' keys + loading_case : str + Loading type + title : str + Plot title (optional) + ax : matplotlib.axes.Axes + Axes to plot on (creates new figure if None) + + Returns + ------- + matplotlib.axes.Axes + The axes with the plot + """ + if ax is None: + fig, ax = plt.subplots(figsize=(8, 6)) + + # Model prediction on fine grid + lambda_model = np.linspace(1.0, lambda_exp.max() * 1.02, 200) + P_model = mooney_rivlin_pk1_stress( + params["C10"], params["C01"], params["kappa"], lambda_model, loading_case + ) + + # Experimental data + ax.plot( + lambda_exp, + P_exp, + "o", + markersize=8, + markerfacecolor="red", + markeredgecolor="black", + label="Treloar experimental", + ) + + # Model prediction + ax.plot( + lambda_model, + P_model, + "-", + linewidth=2, + color="blue", + label=f"Mooney-Rivlin (C10={params['C10']:.4f}, C01={params['C01']:.4f})", + ) + + ax.set_xlabel(r"Stretch $\lambda$ [-]", fontsize=12) + ax.set_ylabel(r"PK1 Stress $P_{11}$ [MPa]", fontsize=12) + ax.set_title(title or f"{loading_case} - Parameter Identification", fontsize=14) + ax.legend(loc="upper left", fontsize=10) + ax.grid(True, alpha=0.3) + + return ax + + +################################################################################### +# Main Execution +# -------------- +# +# We demonstrate two identification strategies: +# +# 1. **Individual fitting**: Optimize parameters separately for each loading case +# 2. **Combined fitting**: Find a single parameter set that best fits all cases +# +# The individual fits provide loading-case-specific parameters (useful for +# understanding material behavior), while combined fitting gives a universal +# parameter set for general use. + +if __name__ == "__main__": + # Locate data file relative to script + # Try multiple approaches to find the data file (supports both direct + # execution and Sphinx-Gallery which does not define __file__) + data_filename = os.path.join("hyperelasticity", "comparison", "Treloar.txt") + possible_paths = [] + + # Try using __file__ if available + try: + script_dir = os.path.dirname(os.path.abspath(__file__)) + possible_paths.append(os.path.join(script_dir, "..", data_filename)) + except NameError: + pass + + # Try relative to current working directory (Sphinx-Gallery context) + possible_paths.append(os.path.join(os.getcwd(), "..", data_filename)) + possible_paths.append(os.path.join(os.getcwd(), "examples", data_filename)) + + # Find first existing path + data_path = None + for path in possible_paths: + if os.path.exists(path): + data_path = path + break + + if data_path is None: + raise FileNotFoundError( + f"Could not find Treloar.txt data file. Searched: {possible_paths}" + ) + + print("=" * 70) + print(" MOONEY-RIVLIN PARAMETER IDENTIFICATION") + print(" Using scipy.differential_evolution and sklearn.mean_squared_error") + print("=" * 70) + + # ------------------------------------------------------------------------- + # Load and prepare experimental data + # ------------------------------------------------------------------------- + df = load_treloar_data(data_path) + print(f"\nLoaded Treloar data: {len(df)} data points") + + # Extract data for each loading case (remove NaN values) + ut_mask = ~df["lambda_1"].isna() & ~df["P1_MPa"].isna() + ps_mask = ~df["lambda_2"].isna() & ~df["P2_MPa"].isna() + et_mask = ~df["lambda_3"].isna() & ~df["P3_MPa"].isna() + + lambda_ut = df.loc[ut_mask, "lambda_1"].values + P_ut = df.loc[ut_mask, "P1_MPa"].values + + lambda_ps = df.loc[ps_mask, "lambda_2"].values + P_ps = df.loc[ps_mask, "P2_MPa"].values + + lambda_et = df.loc[et_mask, "lambda_3"].values + P_et = df.loc[et_mask, "P3_MPa"].values + + print(f"\nData summary:") + print( + f" Uniaxial Tension (UT): {len(lambda_ut)} points, " + f"lambda in [{lambda_ut.min():.2f}, {lambda_ut.max():.2f}]" + ) + print( + f" Pure Shear (PS): {len(lambda_ps)} points, " + f"lambda in [{lambda_ps.min():.2f}, {lambda_ps.max():.2f}]" + ) + print( + f" Equibiaxial (ET): {len(lambda_et)} points, " + f"lambda in [{lambda_et.min():.2f}, {lambda_et.max():.2f}]" + ) + + # Fixed bulk modulus (near-incompressible rubber) + kappa = 4000.0 + + # ------------------------------------------------------------------------- + # Individual fitting for each loading case + # ------------------------------------------------------------------------- + print("\n" + "=" * 70) + print(" INDIVIDUAL FITTING (one parameter set per loading case)") + print("=" * 70) + + params_ut = identify_parameters(lambda_ut, P_ut, kappa, "UT", verbose=True) + params_ps = identify_parameters(lambda_ps, P_ps, kappa, "PS", verbose=True) + params_et = identify_parameters(lambda_et, P_et, kappa, "ET", verbose=True) + + # Reference values from literature (Steinmann et al., 2012) + print("\n" + "-" * 70) + print("Comparison with literature (Steinmann et al., 2012):") + print("-" * 70) + print( + f"{'Case':<8} {'C10 (lit.)':<12} {'C10 (ident.)':<14} {'C01 (lit.)':<12} {'C01 (ident.)':<14}" + ) + print("-" * 70) + print( + f"{'UT':<8} {0.2588:<12.4f} {params_ut['C10']:<14.4f} {-0.0449:<12.4f} {params_ut['C01']:<14.4f}" + ) + print( + f"{'PS':<8} {0.2348:<12.4f} {params_ps['C10']:<14.4f} {-0.065:<12.4f} {params_ps['C01']:<14.4f}" + ) + print( + f"{'ET':<8} {0.1713:<12.4f} {params_et['C10']:<14.4f} {0.0047:<12.4f} {params_et['C01']:<14.4f}" + ) + + # ------------------------------------------------------------------------- + # Combined fitting + # ------------------------------------------------------------------------- + print("\n" + "=" * 70) + print(" COMBINED FITTING (single parameter set for all loading cases)") + print("=" * 70) + + data_combined = { + "UT": (lambda_ut, P_ut), + "PS": (lambda_ps, P_ps), + "ET": (lambda_et, P_et), + } + + params_combined = identify_parameters_combined( + data_combined, kappa=kappa, normalize=True, verbose=True + ) + + # ------------------------------------------------------------------------- + # Visualization + # ------------------------------------------------------------------------- + print("\n" + "=" * 70) + print(" GENERATING PLOTS") + print("=" * 70) + + fig, axes = plt.subplots(2, 3, figsize=(15, 10)) + + # Row 1: Individual fits + plot_fit_comparison( + lambda_ut, P_ut, params_ut, "UT", title="UT - Individual fit", ax=axes[0, 0] + ) + plot_fit_comparison( + lambda_ps, P_ps, params_ps, "PS", title="PS - Individual fit", ax=axes[0, 1] + ) + plot_fit_comparison( + lambda_et, P_et, params_et, "ET", title="ET - Individual fit", ax=axes[0, 2] + ) + + # Row 2: Combined fit applied to all cases + plot_fit_comparison( + lambda_ut, P_ut, params_combined, "UT", title="UT - Combined fit", ax=axes[1, 0] + ) + plot_fit_comparison( + lambda_ps, P_ps, params_combined, "PS", title="PS - Combined fit", ax=axes[1, 1] + ) + plot_fit_comparison( + lambda_et, P_et, params_combined, "ET", title="ET - Combined fit", ax=axes[1, 2] + ) + + fig.suptitle( + "Mooney-Rivlin Parameter Identification using Treloar Data\n" + "(Top: Individual fits, Bottom: Combined fit)", + fontsize=14, + fontweight="bold", + ) + plt.tight_layout() + + # ------------------------------------------------------------------------- + # Summary + # ------------------------------------------------------------------------- + print("\n" + "=" * 70) + print(" SUMMARY OF IDENTIFIED PARAMETERS") + print("=" * 70) + print(f"{'Case':<12} {'C10 [MPa]':>12} {'C01 [MPa]':>12} {'MSE [MPa^2]':>14}") + print("-" * 52) + print( + f"{'UT':<12} {params_ut['C10']:>12.4f} {params_ut['C01']:>12.4f} {params_ut['mse']:>14.2e}" + ) + print( + f"{'PS':<12} {params_ps['C10']:>12.4f} {params_ps['C01']:>12.4f} {params_ps['mse']:>14.2e}" + ) + print( + f"{'ET':<12} {params_et['C10']:>12.4f} {params_et['C01']:>12.4f} {params_et['mse']:>14.2e}" + ) + print("-" * 52) + print( + f"{'Combined':<12} {params_combined['C10']:>12.4f} {params_combined['C01']:>12.4f} {params_combined['cost']:>14.2e}" + ) + print("=" * 70) + + print("\nNote: The combined fit uses NMSE (normalized MSE) which explains") + print("the different cost magnitude compared to individual MSE values.") + + plt.show() diff --git a/include/simcoon/Continuum_mechanics/Homogenization/cylinder_multi.hpp b/include/simcoon/Continuum_mechanics/Homogenization/cylinder_multi.hpp index 444dfd9c8..95f02327d 100755 --- a/include/simcoon/Continuum_mechanics/Homogenization/cylinder_multi.hpp +++ b/include/simcoon/Continuum_mechanics/Homogenization/cylinder_multi.hpp @@ -57,10 +57,6 @@ namespace simcoon{ * - Long fiber reinforced polymers * - Carbon fiber composites * - * **Coordinate System:** - * - The cylinder axis is aligned with the local x₃ direction - * - The cross-section is circular in the x₁-x₂ plane - * * **Concentration Tensors:** * * The strain concentration relates macroscopic to local strain: diff --git a/include/simcoon/Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.hpp b/include/simcoon/Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.hpp index dd4d80798..c9a8fef1d 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.hpp @@ -21,6 +21,7 @@ #pragma once +#include #include namespace simcoon{ @@ -42,7 +43,7 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_hypoelasticity_ortho(const arma::vec &, const arma::vec &, const arma::mat &, const arma::mat&, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &); +void umat_hypoelasticity_ortho(const std::string &umat_name, const arma::vec &etot, const arma::vec &Detot, const arma::mat &F0, const arma::mat &F1, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_finite group diff --git a/include/simcoon/Continuum_mechanics/Umat/Finite/mooney_rivlin.hpp b/include/simcoon/Continuum_mechanics/Umat/Finite/mooney_rivlin.hpp deleted file mode 100755 index 909fd7857..000000000 --- a/include/simcoon/Continuum_mechanics/Umat/Finite/mooney_rivlin.hpp +++ /dev/null @@ -1,50 +0,0 @@ -/* This file is part of simcoon. - - simcoon is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - simcoon is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with simcoon. If not, see . - - */ - -///@file mooney_rivlin.hpp -///@brief User subroutine for Isotropic elastic materials in 3D case -///@version 1.0 - -#pragma once - -#include - -namespace simcoon{ - -/** - * @file mooney_rivlin.hpp - * @brief Finite strain constitutive model. - */ - -/** @addtogroup umat_finite - * @{ - */ - - -///@brief The elastic UMAT requires 2 constants: -///@brief props[0] : Young modulus -///@brief props[1] : Poisson ratio -///@brief props[2] : CTE - -///@brief No statev is required for thermoelastic constitutive law - -void umat_mooney_rivlin(const arma::vec &, const arma::vec &, const arma::mat &, const arma::mat&, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &); - - -/** @} */ // end of umat_finite group - -} //namespace simcoon diff --git a/include/simcoon/Continuum_mechanics/Umat/Finite/neo_hookean_comp.hpp b/include/simcoon/Continuum_mechanics/Umat/Finite/neo_hookean_comp.hpp index 0fe1aaf1d..7ff4f6d39 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Finite/neo_hookean_comp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Finite/neo_hookean_comp.hpp @@ -21,6 +21,7 @@ #pragma once +#include #include namespace simcoon{ @@ -168,7 +169,7 @@ namespace simcoon{ * - Holzapfel, G. A. (2000). *Nonlinear Solid Mechanics: A Continuum Approach for Engineering*. Wiley. * - Connolly, S. J., et al. (2019). "Automatic differentiation based formulation of computational models." *Computational Mechanics*, 64, 1273-1288. */ -void umat_neo_hookean_comp(const arma::vec &Etot, const arma::vec &DEtot, const arma::mat &F0, const arma::mat &F1, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_neo_hookean_comp(const std::string &umat_name, const arma::vec &etot, const arma::vec &Detot, const arma::mat &F0, const arma::mat &F1, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_finite group diff --git a/include/simcoon/Continuum_mechanics/Umat/Finite/neo_hookean_incomp.hpp b/include/simcoon/Continuum_mechanics/Umat/Finite/neo_hookean_incomp.hpp index 03f61b064..c87c11446 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Finite/neo_hookean_incomp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Finite/neo_hookean_incomp.hpp @@ -21,6 +21,7 @@ #pragma once +#include #include namespace simcoon{ @@ -42,7 +43,7 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_neo_hookean_incomp(const arma::vec &etot, const arma::vec &Detot, const arma::mat &F0, const arma::mat &F1, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT,const double &Time,const double &DTime, double &Wm_0, double &Wm_1, double &Wm_2, double &Wm_3, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); +void umat_neo_hookean_incomp(const std::string &umat_name, const arma::vec &etot, const arma::vec &Detot, const arma::mat &F0, const arma::mat &F1, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_finite group diff --git a/include/simcoon/Continuum_mechanics/Umat/Finite/saint_venant.hpp b/include/simcoon/Continuum_mechanics/Umat/Finite/saint_venant.hpp index 50a4d9682..3559bad16 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Finite/saint_venant.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Finite/saint_venant.hpp @@ -21,6 +21,7 @@ #pragma once +#include #include namespace simcoon{ @@ -42,7 +43,7 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_saint_venant(const arma::vec &, const arma::vec &, const arma::mat &, const arma::mat&, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &); +void umat_saint_venant(const std::string &umat_name, const arma::vec &etot, const arma::vec &Detot, const arma::mat &F0, const arma::mat &F1, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_finite group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Damage/damage_LLD_0.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Damage/damage_LLD_0.hpp index 2787fe5c3..45b1e660e 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Damage/damage_LLD_0.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Damage/damage_LLD_0.hpp @@ -24,6 +24,7 @@ #pragma once +#include #include #include @@ -272,7 +273,7 @@ namespace simcoon{ * - Matzenmiller, A., Lubliner, J., & Taylor, R. L. (1995). "A constitutive model for anisotropic damage in fiber-composites." *Mechanics of Materials*, 20(2), 125-152. * - Pinho, S. T., et al. (2012). "Material and structural response of polymer-matrix fibre-reinforced composites." *Journal of Composite Materials*, 46(19-20), 2313-2341. */ -void umat_damage_LLD_0(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_damage_LLD_0(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.hpp index e3f8bb71b..4ebcea6e0 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.hpp @@ -21,6 +21,7 @@ #pragma once +#include #include #include @@ -36,7 +37,7 @@ namespace simcoon{ */ -void umat_damage_weibull(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &); +void umat_damage_weibull(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_isotropic.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_isotropic.hpp index 91d7c6610..957dad9fd 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_isotropic.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_isotropic.hpp @@ -21,6 +21,7 @@ #pragma once +#include #include namespace simcoon{ @@ -122,12 +123,12 @@ namespace simcoon{ * * Total state variables required: \f$ n_{statev} = 1 \f$ * + * @param umat_name Name of the constitutive model (ELISO) * @param Etot Total strain tensor at beginning of increment (Voigt notation: 6×1 vector) * @param DEtot Strain increment tensor (Voigt notation: 6×1 vector) * @param sigma Stress tensor (Voigt notation: 6×1 vector) [output] * @param Lt Tangent modulus \f$ \mathbf{L}_t = \mathbf{L} \f$ (6×6 matrix) [output] * @param L Elastic stiffness tensor (6×6 matrix) [output] - * @param sigma_in Internal stress contribution for explicit solvers (6×1 vector) [output] * @param DR Rotation increment matrix (3×3) for objective integration * @param nprops Number of material properties * @param props Material properties vector (see table above) @@ -144,7 +145,6 @@ namespace simcoon{ * @param ndi Number of direct stress components (typically 3) * @param nshr Number of shear stress components (typically 3) * @param start Flag indicating first increment (true) or continuation (false) - * @param solver_type Solver type: 0=implicit, 1=explicit, 2=dynamic implicit * @param tnew_dt Suggested new time step size for adaptive time stepping [output] * * @note This model is restricted to small strains (< 1%) @@ -169,12 +169,11 @@ namespace simcoon{ * vec sigma = zeros(6); * mat Lt = zeros(6,6); * mat L = zeros(6,6); - * vec sigma_in = zeros(6); * mat DR = eye(3,3); * - * umat_elasticity_iso(Etot, DEtot, sigma, Lt, L, sigma_in, DR, + * umat_elasticity_iso("ELISO", Etot, DEtot, sigma, Lt, L, DR, * 3, props, 1, statev, 25.0, 5.0, 0.0, 1.0, - * Wm, Wm_r, Wm_ir, Wm_d, 3, 3, false, 0, tnew_dt); + * Wm, Wm_r, Wm_ir, Wm_d, 3, 3, false, tnew_dt); * * // Expected stress: sigma(0) ≈ E*Etot(0)*(1-nu)/((1+nu)(1-2nu)) - E*alpha*DT/(1-2nu) * @endcode @@ -184,7 +183,7 @@ namespace simcoon{ * - Bower, A. F. (2009). *Applied Mechanics of Solids*. CRC Press. * - Gurtin, M. E. (1981). *An Introduction to Continuum Mechanics*. Academic Press. */ -void umat_elasticity_iso(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_elasticity_iso(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_orthotropic.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_orthotropic.hpp index 77b817a60..9f38548a6 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_orthotropic.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_orthotropic.hpp @@ -21,6 +21,7 @@ #pragma once +#include #include namespace simcoon{ @@ -97,12 +98,12 @@ namespace simcoon{ * |-------|--------|-------------|-------| * | statev[0] | \f$ T_{init} \f$ | Initial/reference temperature | Temperature | * + * @param umat_name Name of the constitutive model (ELORT) * @param Etot Total strain tensor at beginning of increment (Voigt notation: 6×1) * @param DEtot Strain increment tensor (Voigt notation: 6×1) * @param sigma Stress tensor [output] (Voigt notation: 6×1) * @param Lt Tangent modulus (= L for linear elasticity) [output] (6×6) * @param L Elastic stiffness tensor [output] (6×6) - * @param sigma_in Internal stress for explicit solvers [output] (6×1) * @param DR Rotation increment matrix (3×3) * @param nprops Number of material properties (12) * @param props Material properties vector @@ -119,7 +120,6 @@ namespace simcoon{ * @param ndi Number of direct stress components * @param nshr Number of shear stress components * @param start Flag indicating first increment - * @param solver_type Solver type (0=implicit, 1=explicit) * @param tnew_dt Suggested new time step [output] * * @note Material axes must be aligned with global axes. For rotated materials, @@ -128,7 +128,7 @@ namespace simcoon{ * * @see L_ortho() for orthotropic stiffness tensor construction */ -void umat_elasticity_ortho(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_elasticity_ortho(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_transverse_isotropic.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_transverse_isotropic.hpp index 1f38f2df2..4124510dc 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_transverse_isotropic.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_transverse_isotropic.hpp @@ -21,6 +21,7 @@ #pragma once +#include #include namespace simcoon{ @@ -101,12 +102,12 @@ namespace simcoon{ * |-------|--------|-------------|-------| * | statev[0] | \f$ T_{init} \f$ | Initial/reference temperature | Temperature | * + * @param umat_name Name of the constitutive model (ELIST) * @param Etot Total strain tensor at beginning of increment (Voigt notation: 6×1) * @param DEtot Strain increment tensor (Voigt notation: 6×1) * @param sigma Stress tensor [output] (Voigt notation: 6×1) * @param Lt Tangent modulus (= L for linear elasticity) [output] (6×6) * @param L Elastic stiffness tensor [output] (6×6) - * @param sigma_in Internal stress for explicit solvers [output] (6×1) * @param DR Rotation increment matrix (3×3) * @param nprops Number of material properties (7) * @param props Material properties vector @@ -123,7 +124,6 @@ namespace simcoon{ * @param ndi Number of direct stress components * @param nshr Number of shear stress components * @param start Flag indicating first increment - * @param solver_type Solver type (0=implicit, 1=explicit) * @param tnew_dt Suggested new time step [output] * * @note The fiber direction (axis 1) must be aligned with the global x-axis. @@ -132,7 +132,7 @@ namespace simcoon{ * * @see L_isotrans() for transversely isotropic stiffness tensor construction */ -void umat_elasticity_trans_iso(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_elasticity_trans_iso(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/External/external_umat.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/External/external_umat.hpp index 2eb10790d..9bab02edd 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/External/external_umat.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/External/external_umat.hpp @@ -1,26 +1,176 @@ /* This file is part of simcoon. - + simcoon is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. - + simcoon is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - + You should have received a copy of the GNU General Public License along with simcoon. If not, see . - + */ ///@file external_umat.hpp -///@brief User subroutine for Isotropic elastic materials in 3D case +///@brief External user material subroutine interface for custom mechanical models ///@version 1.0 #pragma once #include -extern void umat_external(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &); +/** + * @file external_umat.hpp + * @brief External user material subroutine interface for custom mechanical models + * @author Yves Chemisky + * @version 1.0 + * + * This file declares the external interface for user-defined mechanical constitutive + * models that can be loaded as plugins at runtime. Users implement their material + * models through the plugin API system. + */ + +/** @addtogroup umat_mechanical + * @{ + */ + +/** + * @brief External mechanical UMAT subroutine interface + * + * @details This function declaration provides the interface for external user-defined + * material subroutines (UMATs). The actual implementation is provided by the user + * through the plugin system, allowing custom constitutive models to be loaded + * dynamically at runtime. + * + * **Plugin System:** + * + * To implement a custom material model, users should: + * 1. Create a class inheriting from `umat_plugin_ext_api` + * 2. Implement the `umat_external_M()` method with the constitutive equations + * 3. Compile as a shared library (`.so`, `.dylib`, or `.dll`) + * 4. Configure the material name in `material.dat` to match the plugin's `name()` + * + * **Available Plugin Formats:** + * + * | Format | Class | Description | + * |--------|-------|-------------| + * | UMEXT | `umat_plugin_ext_api` | Native simcoon format using Armadillo | + * | UMABA | `umat_plugin_aba_api` | Abaqus UMAT-compatible format | + * | UMANS | `umat_plugin_ans_api` | ANSYS USERMAT-compatible format | + * + * **Implementation Example:** + * + * @code + * #include + * + * class MyMaterial : public umat_plugin_ext_api { + * public: + * std::string name() const override { return "my_material"; } + * + * void umat_external_M( + * const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, + * arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, + * const arma::mat &DR, const int &nprops, + * const arma::vec &props, const int &nstatev, arma::vec &statev, + * const double &T, const double &DT, const double &Time, const double &DTime, + * double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, + * const int &ndi, const int &nshr, const bool &start, + * const int &solver_type, double &tnew_dt + * ) override { + * // Implement your constitutive model here + * // 1. Extract material properties from props + * // 2. Update state variables in statev + * // 3. Compute stress tensor sigma + * // 4. Compute tangent modulus Lt (for implicit solvers) + * // 5. Update work quantities Wm, Wm_r, Wm_ir, Wm_d + * } + * }; + * + * // Required factory functions for plugin loading + * extern "C" umat_plugin_ext_api* create_api() { + * return new MyMaterial(); + * } + * + * extern "C" void destroy_api(umat_plugin_ext_api* p) { + * delete p; + * } + * @endcode + * + * **Compiling the Plugin:** + * + * @code{.sh} + * # On macOS + * clang++ -c -fPIC -std=c++14 my_material.cpp -I/path/to/simcoon/include + * clang++ -std=c++14 -shared -lsimcoon -larmadillo -o libmy_material.dylib my_material.o + * + * # On Linux + * g++ -c -fPIC -std=c++14 my_material.cpp -I/path/to/simcoon/include + * g++ -std=c++14 -shared -lsimcoon -larmadillo -o libmy_material.so my_material.o + * @endcode + * + * **Material Configuration (material.dat):** + * + * @code{.txt} + * Material + * Name MY_MATERIAL + * Number_of_material_parameters 3 + * Number_of_internal_variables 1 + * + * #Mechanical + * E 210000 + * nu 0.3 + * alpha 1.2e-5 + * @endcode + * + * **Material Parameters (props):** + * + * User-defined. The number and meaning of material parameters depends on + * the specific constitutive model implemented in the plugin. + * + * **State Variables (statev):** + * + * User-defined. The number and meaning of state variables depends on + * the specific constitutive model implemented in the plugin. + * + * @param Etot Total strain tensor at beginning of increment (Voigt notation: 6×1 vector) + * @param DEtot Strain increment tensor (Voigt notation: 6×1 vector) + * @param sigma Stress tensor (Voigt notation: 6×1 vector) [output] + * @param Lt Consistent tangent modulus (6×6 matrix) [output] + * @param L Elastic stiffness tensor (6×6 matrix) [output] + * @param sigma_in Internal stress contribution for explicit solvers (6×1 vector) [output] + * @param DR Rotation increment matrix (3×3) for objective integration + * @param nprops Number of material properties + * @param props Material properties vector + * @param nstatev Number of state variables + * @param statev State variables vector [input/output] + * @param T Temperature at beginning of increment + * @param DT Temperature increment + * @param Time Time at beginning of increment + * @param DTime Time increment + * @param Wm Total mechanical work [output] + * @param Wm_r Recoverable (elastic) work [output] + * @param Wm_ir Irrecoverable work [output] + * @param Wm_d Dissipated work [output] + * @param ndi Number of direct stress components (typically 3) + * @param nshr Number of shear stress components (typically 3) + * @param start Flag indicating first increment (true) or continuation (false) + * @param solver_type Solver type: 0=implicit (Newton), 1=explicit (RNL), 2=dynamic implicit + * @param tnew_dt Suggested new time step ratio for adaptive time stepping [output] + * + * @note This is an external declaration; the function must be provided by a loaded plugin + * @note The plugin shared library must be in the runtime library path + * @note Use `start` flag to initialize state variables on first call + * @note For implicit solvers (solver_type=0,2), compute consistent tangent Lt + * @note For explicit solvers (solver_type=1), compute internal stress sigma_in + * + * @see umat_plugin_ext_api for the plugin base class interface + * @see umat_plugin_aba_api for Abaqus-compatible plugins + * @see umat_plugin_ans_api for ANSYS-compatible plugins + */ +extern void umat_external(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); + +/** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.hpp index b88834527..1211dd36e 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon{ @@ -142,7 +143,7 @@ namespace simcoon{ * - Chaboche, J. L. (1986). "Time-independent constitutive theories for cyclic plasticity." *Int. J. Plasticity*, 2(2), 149-188. * - Barlat, F., et al. (2005). "Linear transformation-based anisotropic yield functions." *Int. J. Plasticity*, 21(5), 1009-1039. */ -void umat_ani_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_ani_chaboche_CCP(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.hpp index ff2ff3af1..3712c7d8d 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon{ @@ -129,7 +130,7 @@ namespace simcoon{ * - Teodosiu, C., & Hu, Z. (1998). "Microstructure in the continuum modeling of plastic anisotropy." *Proc. 19th Riso Int. Symp.*, 149-168. * - Holmedal, B. (2019). "Bauschinger effect modelled by yield surface distortions." *Int. J. Plasticity*, 123, 86-100. */ -void umat_dfa_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_dfa_chaboche_CCP(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.hpp index 0d365d48d..7c910b4a7 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon{ @@ -103,7 +104,7 @@ namespace simcoon{ * @param sigma Stress tensor (Voigt notation: 6×1 vector) [output] * @param Lt Consistent tangent modulus (6×6 matrix) [output] * @param L Elastic stiffness tensor (6×6 matrix) [output] - * @param sigma_in Internal stress for explicit solvers (6×1 vector) [output] + * @param umat_name Name of the constitutive model * @param DR Rotation increment matrix (3×3) for objective integration * @param nprops Number of material properties * @param props Material properties vector (criterion-dependent) @@ -120,7 +121,6 @@ namespace simcoon{ * @param ndi Number of direct stress components (typically 3) * @param nshr Number of shear stress components (typically 3) * @param start Flag indicating first increment (true) or continuation (false) - * @param solver_type Solver type: 0=implicit, 1=explicit, 2=dynamic implicit * @param tnew_dt Suggested new time step size for adaptive time stepping [output] * * @note Use this as a base for implementing new yield criteria @@ -135,7 +135,7 @@ namespace simcoon{ * - Chaboche, J. L. (2008). "A review of some plasticity and viscoplasticity constitutive theories." *Int. J. Plasticity*, 24(10), 1642-1693. * - Simo, J. C., & Hughes, T. J. R. (1998). *Computational Inelasticity*. Springer. */ -void umat_generic_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_generic_chaboche_CCP(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.hpp index e64b8c555..01b2e4456 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon{ @@ -154,7 +155,7 @@ namespace simcoon{ * - Chaboche, J. L. (1986). "Time-independent constitutive theories for cyclic plasticity." *Int. J. Plasticity*, 2(2), 149-188. * - Barlat, F., et al. (2005). "Linear transformation-based anisotropic yield functions." *Int. J. Plasticity*, 21(5), 1009-1039. */ -void umat_hill_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_hill_chaboche_CCP(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.hpp index 1640c6d14..b0efadf07 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon { @@ -257,7 +258,7 @@ namespace simcoon { * - Banabic, D. (2010). *Sheet Metal Forming Processes: Constitutive Modelling and Numerical Simulation*. Springer. * - Ortiz, M., & Simo, J. C. (1986). "An analysis of a new class of integration algorithms for elastoplastic constitutive relations." *International Journal for Numerical Methods in Engineering*, 23(3), 353-366. */ -void umat_plasticity_hill_isoh_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); +void umat_plasticity_hill_isoh_CCP(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.hpp index c2174ae53..7bc825ffc 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon { @@ -147,7 +148,7 @@ namespace simcoon { * - Voce, E. (1948). "The relationship between stress and strain for homogeneous deformation." *J. Inst. Met.*, 74, 537-562. * - Barlat, F., et al. (2003). "Plane stress yield function for aluminum alloy sheets." *Int. J. Plasticity*, 19(9), 1297-1319. */ -void umat_plasticity_hill_isoh_CCP_N(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); +void umat_plasticity_hill_isoh_CCP_N(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.hpp index a85e77e7f..16ebd2f88 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon{ @@ -262,7 +263,7 @@ namespace simcoon{ * - Lemaitre, J., & Chaboche, J. L. (1990). *Mechanics of Solid Materials*. Cambridge University Press. * - Ortiz, M., & Simo, J. C. (1986). "An analysis of a new class of integration algorithms for elastoplastic constitutive relations." *International Journal for Numerical Methods in Engineering*, 23(3), 353-366. */ -void umat_plasticity_chaboche_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_plasticity_chaboche_CCP(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.hpp index a9966689c..364009cac 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.hpp @@ -22,6 +22,7 @@ ///@version 1.0 #pragma once +#include #include namespace simcoon{ @@ -104,12 +105,12 @@ namespace simcoon{ * | statev[6] | \f$ \varepsilon^p_{13} \f$ | Plastic strain component 13 (×2 in Voigt) | - | * | statev[7] | \f$ \varepsilon^p_{23} \f$ | Plastic strain component 23 (×2 in Voigt) | - | * + * @param umat_name Name of the constitutive model (EPICP) * @param Etot Total strain tensor at beginning of increment (Voigt notation: 6×1 vector) * @param DEtot Strain increment tensor (Voigt notation: 6×1 vector) * @param sigma Stress tensor (Voigt notation: 6×1 vector) [output] * @param Lt Consistent tangent modulus \f$ \mathbf{L}_t = \frac{\partial \boldsymbol{\sigma}}{\partial \boldsymbol{\varepsilon}} \f$ (6×6 matrix) [output] * @param L Elastic stiffness tensor (6×6 matrix) [output] - * @param sigma_in Internal stress contribution for explicit solvers (6×1 vector) [output] * @param DR Rotation increment matrix (3×3) for objective integration * @param nprops Number of material properties * @param props Material properties vector (see table above) @@ -126,12 +127,11 @@ namespace simcoon{ * @param ndi Number of direct stress components (typically 3) * @param nshr Number of shear stress components (typically 3) * @param start Flag indicating first increment (true) or continuation (false) - * @param solver_type Solver type: 0=implicit, 1=explicit, 2=dynamic implicit * @param tnew_dt Suggested new time step size for adaptive time stepping [output] * * @note Voigt notation convention: [11, 22, 33, 12, 13, 23] with engineering shear strains (γ = 2ε) * @note The consistent tangent modulus Lt ensures quadratic convergence in implicit Newton-Raphson schemes - * @note For explicit solvers (solver_type=1), use sigma_in instead of Lt + * @note The tangent modulus Lt is always computed * * @see Fischer_Burmeister_m() for the complementarity solver * @see L_iso() for isotropic elastic stiffness construction @@ -144,17 +144,16 @@ namespace simcoon{ * vec sigma = zeros(6); * mat Lt = zeros(6,6); * mat L = zeros(6,6); - * vec sigma_in = zeros(6); * mat DR = eye(3,3); * vec props = {70000, 0.3, 1e-5, 200, 500, 0.2}; * vec statev = zeros(8); * - * umat_plasticity_iso_CCP(Etot, DEtot, sigma, Lt, L, sigma_in, DR, 6, props, 8, statev, + * umat_plasticity_iso_CCP("EPICP", Etot, DEtot, sigma, Lt, L, DR, 6, props, 8, statev, * 20.0, 0.0, 0.0, 1.0, Wm, Wm_r, Wm_ir, Wm_d, - * 3, 3, true, 0, tnew_dt); + * 3, 3, true, tnew_dt); * @endcode */ -void umat_plasticity_iso_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_plasticity_iso_CCP(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.hpp index 42d29bf6a..65756700c 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon{ @@ -169,7 +170,7 @@ namespace simcoon{ * - Lemaitre, J., & Chaboche, J. L. (1990). *Mechanics of Solid Materials*. Cambridge University Press. * - Simo, J. C., & Hughes, T. J. R. (1998). *Computational Inelasticity*. Springer. */ -void umat_plasticity_kin_iso_CCP(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt); +void umat_plasticity_kin_iso_CCP(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.hpp index 7579b7ca8..139d20391 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.hpp @@ -23,6 +23,7 @@ */ #pragma once +#include #include namespace simcoon{ @@ -39,12 +40,24 @@ namespace simcoon{ * crystal by tracking N individual martensite variants, where N is typically 12 or 24 depending * on the crystallographic system. * + * **Elastic Symmetry Selection:** + * + * The elastic behavior is determined by the umat_name parameter: + * + * | umat_name | Symmetry | Elastic Parameters | + * |-----------|----------|-------------------| + * | SMAMO | Isotropic | E, nu | + * | SMAMC | Cubic | C11, C12, C44 | + * | SMAOT | Orthotropic | E1, E2, E3, nu12, nu13, nu23, G12, G13, G23 | + * | SMATI | Transverse Isotropic | EL, ET, nuTL, nuTT, GLT | + * * **Key Features:** * - Explicit tracking of N martensite variant volume fractions \f$ f_n \f$ (n = 1, ..., N) * - Crystallographic transformation strains for each variant * - Variant selection driven by resolved thermodynamic driving force * - Inter-variant interactions through hardening matrix * - Applicable to single crystal SMA behavior + * - Supports multiple elastic symmetry types * * **Physical Background:** * @@ -121,9 +134,40 @@ namespace simcoon{ * - Cubic → Monoclinic (NiTi): N = 12 or 24 variants * - Cubic → Tetragonal: N = 3 variants * + * **Material Parameters (props):** + * + * The props vector layout depends on the elastic symmetry (umat_name): + * + * **SMAMO (Isotropic):** 16 parameters + * | Index | Symbol | Description | + * |-------|--------|-------------| + * | props[0] | E | Young's modulus | + * | props[1] | nu | Poisson's ratio | + * | props[2] | alpha_iso | Coefficient of thermal expansion | + * | props[3] | b | Slope parameter | + * | props[4] | g | Shear strain magnitude | + * | props[5] | Ms | Martensite start temperature | + * | props[6] | Af | Austenite finish temperature | + * | props[7] | nvariants | Number of martensite variants | + * | props[8-15] | c_lambda0, p0_lambda0, ... | Lagrange multiplier parameters | + * + * **SMAMC (Cubic):** 17 parameters + * | Index | Symbol | Description | + * |-------|--------|-------------| + * | props[0] | C11 | Elastic constant C11 | + * | props[1] | C12 | Elastic constant C12 | + * | props[2] | C44 | Elastic constant C44 | + * | props[3] | alpha_iso | Coefficient of thermal expansion | + * | props[4] | b | Slope parameter | + * | props[5] | g | Shear strain magnitude | + * | props[6] | Ms | Martensite start temperature | + * | props[7] | Af | Austenite finish temperature | + * | props[8] | nvariants | Number of martensite variants | + * | props[9-16] | c_lambda0, p0_lambda0, ... | Lagrange multiplier parameters | + * * **State Variables (statev):** * - * Total state variables required: \f$ n_{statev} = 1 + N + 6 \f$ (for N variants) + * Total state variables required: \f$ n_{statev} = 1 + N + 7 \f$ (for N variants) * * | Index | Symbol | Description | Units | * |-------|--------|-------------|-------| @@ -132,16 +176,13 @@ namespace simcoon{ * | statev[2] | \f$ f_2 \f$ | Volume fraction of martensite variant 2 | - | * | ... | ... | ... | ... | * | statev[N] | \f$ f_N \f$ | Volume fraction of martensite variant N | - | - * | statev[N+1] | \f$ \varepsilon^{tr}_{11} \f$ | Macroscopic transformation strain component 11 | Strain | - * | statev[N+2] | \f$ \varepsilon^{tr}_{22} \f$ | Macroscopic transformation strain component 22 | Strain | - * | statev[N+3] | \f$ \varepsilon^{tr}_{33} \f$ | Macroscopic transformation strain component 33 | Strain | - * | statev[N+4] | \f$ \varepsilon^{tr}_{12} \f$ | Macroscopic transformation strain component 12 | Strain | - * | statev[N+5] | \f$ \varepsilon^{tr}_{13} \f$ | Macroscopic transformation strain component 13 | Strain | - * | statev[N+6] | \f$ \varepsilon^{tr}_{23} \f$ | Macroscopic transformation strain component 23 | Strain | + * | statev[N+1:N+6] | \f$ \varepsilon^{tr} \f$ | Macroscopic transformation strain (Voigt) | Strain | + * | statev[N+7] | \f$ \xi \f$ | Total martensite volume fraction | - | * * The macroscopic transformation strain is computed as: * \f$ \boldsymbol{\varepsilon}^{tr} = \sum_{n=1}^{N} f_n \boldsymbol{\varepsilon}^{tr}_n \f$ * + * @param umat_name Name of the constitutive model (SMAMO, SMAMC, SMAOT, SMATI) * @param Etot Total strain tensor at beginning of increment (Voigt notation: 6×1 vector) * @param DEtot Strain increment tensor (Voigt notation: 6×1 vector) * @param sigma Cauchy stress tensor (Voigt notation: 6×1 vector) [output] @@ -149,7 +190,7 @@ namespace simcoon{ * @param L Elastic stiffness tensor (6×6 matrix) [output] * @param DR Rotation increment matrix (3×3) for objective integration * @param nprops Number of material properties - * @param props Material properties vector + * @param props Material properties vector (layout depends on umat_name) * @param nstatev Number of state variables (includes N variant volume fractions) * @param statev State variables vector containing variant fractions [input/output] * @param T Temperature at beginning of increment @@ -170,6 +211,11 @@ namespace simcoon{ * @note The number of variants N depends on the crystallographic transformation system * @note Variant transformation strains must be provided based on crystallographic data * + * @see L_iso() for isotropic stiffness tensor (SMAMO) + * @see L_cubic() for cubic stiffness tensor (SMAMC) + * @see L_ortho() for orthotropic stiffness tensor (SMAOT) + * @see L_isotrans() for transverse isotropic stiffness tensor (SMATI) + * * **References:** * - Patoor, E., Eberhardt, A., & Berveiller, M. (1996). "Micromechanical modelling of * superelasticity in shape memory alloys." *Journal de Physique IV*, 6(C1), 277-292. @@ -179,7 +225,7 @@ namespace simcoon{ * - Gall, K., & Sehitoglu, H. (1999). "The role of texture in tension-compression * asymmetry in polycrystalline NiTi." *International Journal of Plasticity*, 15(1), 69-92. */ -void umat_sma_mono(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); +void umat_sma_mono(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.hpp deleted file mode 100755 index 737a3394b..000000000 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.hpp +++ /dev/null @@ -1,79 +0,0 @@ -/* This file is part of simcoon. - - simcoon is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - simcoon is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with simcoon. If not, see . - - */ - -///@file SMA_mono_cubic.hpp -///@brief Unified model from: -///@brief Constitutive model of Patoor - 1996, implemented by Y. Chemisky and D. Chatziathanasiou - -#pragma once -#include - -namespace simcoon{ - -///@brief The unified SMA UMAT for transformation requires 25 constants and 17 statev: -///@brief The mechanical transformation UMAT for SMAs has the following statev and props - -///@brief props[0] : flagT: 0 transformation temperatures linearly extrapolated; 1 : smooth temperatures -///@brief props[1] : EA: Young's modulus of Austenite -///@brief props[2] : EM: Young's modulus of Martensite -///@brief props[3] : nuA : Poisson's ratio of Austenite -///@brief props[4] : nuM : Poisson's ratio of Martensite -///@brief props[5] : alphaA_iso : CTE of Austenite -///@brief props[6] : alphaM_iso : CTE of Martensite -///@brief props[7] : Hmin : Minimal transformation strain magnitude -///@brief props[8] : Hmax : Maximal transformation strain magnitude -///@brief props[9] : k1 : Exponential evolution of transformation strain magnitude -///@brief props[10] : sigmacrit : Critical stress for change of transformation strain magnitude -///@brief props[11]: C_A : Slope of martesnite -> austenite parameter -///@brief props[12]: C_M : Slope of austenite -> martensite parameter -///@brief props[13]: Ms0 : Martensite start at zero stress -///@brief props[14]: Mf0 : Martensite finish at zero stress -///@brief props[15]: As0 : Austenite start at zero stress -///@brief props[16]: Af0 : Austenite finish at zero stress -///@brief props[17]: n1 : Martensite start smooth exponent -///@brief props[18]: n2 : Martensite finish smooth exponent -///@brief props[19]: n3 : Austenite start smooth exponent -///@brief props[20]: n4 : Austenite finish smooth exponent - -///@brief props[21]: c_lambda : penalty function exponent start point -///@brief props[22]: p0_lambda : penalty function exponent limit penalty value -///@brief props[23]: n_lambda : penalty function power law exponent -///@brief props[24]: alpha_lambda : penalty function power law parameter - -///@brief The elastic-plastic UMAT with isotropic hardening requires 14 statev: -///@brief statev[0] : T_init : Initial temperature -///@brief statev[1] : xi : MVF (Martensitic volume fraction) -///@brief statev[2] : Transformation strain 11: ET(0,0) -///@brief statev[3] : Transformation strain 22: ET(1,1) -///@brief statev[4] : Transformation strain 33: ET(2,2) -///@brief statev[5] : Transformation strain 12: ET(0,1) (*2) -///@brief statev[6] : Transformation strain 13: ET(0,2) (*2) -///@brief statev[7] : Transformation strain 23: ET(1,2) (*2) - -///@brief statev[8] : xiF : forward MVF -///@brief statev[9] : xiR : reverse MVF -///@brief statev[10] : rhoDs0 difference in entropy for the phases (M - A) -///@brief statev[11] : rhoDs0 difference in internal energy for the phases (M - A) -///@brief statev[12] : parameter for the stress dependance of transformation limits -///@brief statev[13] : a1 : forward hardening parameter -///@brief statev[14] : a2 : reverse hardening parameter -///@brief statev[15] : a3 : Equilibrium hardening parameter -///@brief statev[16] : Y0t : Initial transformation critical value - -void umat_sma_mono_cubic(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &); - -} //namespace simcoon diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/aniso_T.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/aniso_T.hpp index 179a05362..8bd1f270a 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/aniso_T.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/aniso_T.hpp @@ -21,6 +21,7 @@ ///@brief Implemented in 1D-2D-3D #pragma once +#include #include namespace simcoon{ @@ -86,6 +87,6 @@ namespace simcoon{ ///@brief statev[15] : a3 : Equilibrium hardening parameter ///@brief statev[16] : Y0t : Initial transformation critical value -void umat_sma_aniso_T(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &); +void umat_sma_aniso_T(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); } //namespace simcoon diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.hpp index 651f68aa5..01744bcff 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.hpp @@ -21,6 +21,7 @@ ///@brief Implemented in 1D-2D-3D #pragma once +#include #include namespace simcoon{ @@ -77,6 +78,6 @@ namespace simcoon{ ///@brief statev[15] : a3 : Equilibrium hardening parameter ///@brief statev[16] : Y0t : Initial transformation critical value -void umat_sma_unified_T(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &); +void umat_sma_unified_T(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); } //namespace simcoon diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.hpp index ac96f9d69..2f4d26eef 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.hpp @@ -24,6 +24,7 @@ #pragma once +#include #include #include @@ -167,7 +168,7 @@ namespace simcoon { * - Simo, J. C., & Hughes, T. J. R. (1998). *Computational Inelasticity*. Springer. * - Park, S. W., & Schapery, R. A. (1999). "Methods of interconversion between linear viscoelastic material functions. Part I—A numerical method based on Prony series." *International Journal of Solids and Structures*, 36(11), 1653-1675. */ -void umat_prony_Nfast(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); +void umat_prony_Nfast(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_Nfast.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_Nfast.hpp index 69ea29ac0..a406668fc 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_Nfast.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_Nfast.hpp @@ -24,6 +24,7 @@ #pragma once +#include #include #include @@ -171,7 +172,7 @@ namespace simcoon { * - Simo, J. C., & Hughes, T. J. R. (1998). *Computational Inelasticity*. Springer. * - Park, S. W., & Schapery, R. A. (1999). "Methods of interconversion between linear viscoelastic material functions." *Int. J. Solids Struct.*, 36(11), 1653-1675. */ -void umat_zener_Nfast(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); +void umat_zener_Nfast(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_fast.hpp b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_fast.hpp index 70fbc0e1b..292c9cfe8 100755 --- a/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_fast.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_fast.hpp @@ -24,6 +24,7 @@ #pragma once +#include #include #include @@ -147,7 +148,7 @@ namespace simcoon { * - Ferry, J. D. (1980). *Viscoelastic Properties of Polymers* (3rd ed.). Wiley. * - Simo, J. C., & Hughes, T. J. R. (1998). *Computational Inelasticity*. Springer. */ -void umat_zener_fast(const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); +void umat_zener_fast(const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt); /** @} */ // end of umat_mechanical group diff --git a/include/simcoon/Continuum_mechanics/Umat/umat_plugin_api.hpp b/include/simcoon/Continuum_mechanics/Umat/umat_plugin_api.hpp index 404486b34..e76e56949 100755 --- a/include/simcoon/Continuum_mechanics/Umat/umat_plugin_api.hpp +++ b/include/simcoon/Continuum_mechanics/Umat/umat_plugin_api.hpp @@ -93,12 +93,12 @@ class umat_plugin_ext_api { * This function is called at each integration point to update stress and * compute the consistent tangent modulus. * + * @param umat_name Name of the constitutive model * @param Etot Total strain tensor at beginning of increment (Voigt notation: 6×1) * @param DEtot Strain increment tensor (Voigt notation: 6×1) * @param sigma Stress tensor [output] (Voigt notation: 6×1) * @param Lt Consistent tangent modulus [output] (6×6) * @param L Elastic stiffness tensor [output] (6×6) - * @param sigma_in Internal stress for explicit solvers [output] (6×1) * @param DR Rotation increment matrix (3×3) * @param nprops Number of material properties * @param props Material properties vector @@ -115,17 +115,16 @@ class umat_plugin_ext_api { * @param ndi Number of direct stress components * @param nshr Number of shear stress components * @param start Flag indicating first increment - * @param solver_type Solver type (0=implicit, 1=explicit, 2=dynamic) * @param tnew_dt Suggested new time step [output] */ virtual void umat_external_M( - const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, - arma::mat &Lt, arma::mat &L, arma::vec &sigma_in, + const std::string &umat_name, const arma::vec &Etot, const arma::vec &DEtot, arma::vec &sigma, + arma::mat &Lt, arma::mat &L, const arma::mat &DR, const int &nprops, const arma::vec &props, const int &nstatev, arma::vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, - const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt + const int &ndi, const int &nshr, const bool &start, double &tnew_dt ) = 0; /** @brief Virtual destructor for proper cleanup */ diff --git a/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/umat.cpp b/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/umat.cpp index 8641b994e..b8ed5c1a0 100644 --- a/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/umat.cpp +++ b/simcoon-python-builder/src/python_wrappers/Libraries/Continuum_mechanics/umat.cpp @@ -29,7 +29,6 @@ #include #include #include -#include #include #include #include @@ -73,13 +72,10 @@ namespace simpy { int id_umat = list_umat[umat_name_py]; int arguments_type; //depends on the argument used in the umat - void (*umat_function)(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &); - void (*umat_function_2)(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &); - void (*umat_function_3)(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &); - void (*umat_function_4)(const std::string &, const arma::vec &, const arma::vec &, const arma::mat &, const arma::mat &, arma::vec &, arma::mat &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &); - - //scalar needed to launch umat - const int solver_type = 0; + // Unified small-strain function pointer: (umat_name, Etot, DEtot, sigma, Lt, L, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt) + void (*umat_function)(const std::string &, const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &, const double &, const double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &); + // Unified finite-strain function pointer: (umat_name, etot, Detot, F0, F1, sigma, Lt, L, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt) + void (*umat_function_finite)(const std::string &, const arma::vec &, const arma::vec &, const arma::mat &, const arma::mat &, arma::vec &, arma::mat &, arma::mat &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &, const double &, const double &, double &, double &, double &, double &, const int &, const int &, const bool &, double &); const int ncomp=6; int nshr; if (ndi==3) { @@ -134,7 +130,6 @@ namespace simpy { mat list_Wm = carma::arr_to_mat(std::move(Wm_py)); //copy data because values are changed by the umat and returned to python cube L(ncomp, ncomp, nb_points); cube Lt(ncomp, ncomp, nb_points); - vec sigma_in = zeros(1); //not used int nprops = list_props.n_rows; int nstatev = list_statev.n_rows; @@ -142,128 +137,102 @@ namespace simpy { case 2: { umat_function = &simcoon::umat_elasticity_iso; arguments_type = 1; - //simcoon::umat_elasticity_iso(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; } case 3: { umat_function = &simcoon::umat_elasticity_trans_iso; arguments_type = 1; - //simcoon::umat_elasticity_trans_iso(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; } case 4: { umat_function = &simcoon::umat_elasticity_ortho; arguments_type = 1; - //simcoon::umat_elasticity_ortho(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; } case 5: { umat_function = &simcoon::umat_plasticity_iso_CCP; arguments_type = 1; - //simcoon::umat_plasticity_iso_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; } case 6: { umat_function = &simcoon::umat_plasticity_kin_iso_CCP; arguments_type = 1; - //simcoon::umat_plasticity_kin_iso_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; } case 7: { umat_function = &simcoon::umat_plasticity_chaboche_CCP; arguments_type = 1; - //simcoon::umat_plasticity_chaboche_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; } case 8: { - umat_function_2 = &simcoon::umat_plasticity_hill_isoh_CCP; - arguments_type = 2; - //simcoon::umat_plasticity_hill_isoh_CCP(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt); + umat_function = &simcoon::umat_plasticity_hill_isoh_CCP; + arguments_type = 1; break; } case 9: { umat_function = &simcoon::umat_hill_chaboche_CCP; arguments_type = 1; - //simcoon::umat_hill_chaboche_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; } case 10: { umat_function = &simcoon::umat_ani_chaboche_CCP; arguments_type = 1; - //simcoon::umat_ani_chaboche_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; } case 11: { umat_function = &simcoon::umat_dfa_chaboche_CCP; arguments_type = 1; - //simcoon::umat_dfa_chaboche_CCP(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; - } + } case 12: { - umat_function_2 = &simcoon::umat_plasticity_hill_isoh_CCP_N; - arguments_type = 2; - //simcoon::umat_plasticity_hill_isoh_CCP_N(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); + umat_function = &simcoon::umat_plasticity_hill_isoh_CCP_N; + arguments_type = 1; break; - } - case 13: { - umat_function_2 = &simcoon::umat_sma_unified_T; - arguments_type = 2; - //simcoon::umat_sma_unified_T(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt); + } + case 13: { + umat_function = &simcoon::umat_sma_unified_T; + arguments_type = 1; break; } - case 14: { - umat_function_2 = &simcoon::umat_sma_aniso_T; - arguments_type = 2; - //simcoon::umat_sma_aniso_T(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt); + case 14: { + umat_function = &simcoon::umat_sma_aniso_T; + arguments_type = 1; break; } case 15: { umat_function = &simcoon::umat_damage_LLD_0; arguments_type = 1; - //simcoon::umat_damage_LLD_0(etot, Detot, sigma, Lt, L, sigma_in, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, solver_type, tnew_dt); break; - } + } case 16: { - umat_function_2 = &simcoon::umat_zener_fast; - arguments_type = 2; - //simcoon::umat_zener_fast(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt); + umat_function = &simcoon::umat_zener_fast; + arguments_type = 1; break; - } + } case 17: { - umat_function_2 = &simcoon::umat_zener_Nfast; - arguments_type = 2; - //simcoon::umat_zener_Nfast(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt); + umat_function = &simcoon::umat_zener_Nfast; + arguments_type = 1; break; } case 18: { - umat_function_2 = &simcoon::umat_prony_Nfast; - arguments_type = 2; - //simcoon::umat_prony_Nfast(etot, Detot, sigma, Lt, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt); - break; - } - case 19: { - umat_function_3 = &simcoon::umat_sma_mono; - arguments_type = 3; - //simcoon::umat_sma_mono(etot, Detot, sigma, Lt, L, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt); + umat_function = &simcoon::umat_prony_Nfast; + arguments_type = 1; break; } - case 20: { - umat_function_3 = &simcoon::umat_sma_mono_cubic; - arguments_type = 3; - //simcoon::umat_sma_mono_cubic(umat_name, etot, Detot, F0, F1, sigma, Lt, L, DR, nprops, props, nstatev, statev, T, DT, Time, DTime, Wm, Wm_r, Wm_ir, Wm_d, ndi, nshr, start, tnew_dt); + case 19: case 20: { + umat_function = &simcoon::umat_sma_mono; + arguments_type = 1; break; } case 21: case 22: case 23: case 24: case 26: { F0 = carma::arr_to_cube_view(F0_py); - F1 = carma::arr_to_cube_view(F1_py); - umat_function_4 = &simcoon::umat_generic_hyper_invariants; - arguments_type = 4; - break; - } + F1 = carma::arr_to_cube_view(F1_py); + umat_function_finite = &simcoon::umat_generic_hyper_invariants; + arguments_type = 2; + break; + } default: { - //py::print("Error: The choice of Umat could not be found in the umat library \n"); throw std::invalid_argument( "The choice of Umat could not be found in the umat library." ); - //exit(0); } } @@ -294,21 +263,13 @@ namespace simpy { switch (arguments_type) { case 1: { - umat_function(etot, Detot, sigma, Lt.slice(pt), L.slice(pt), sigma_in, DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_function(umat_name_py, etot, Detot, sigma, Lt.slice(pt), L.slice(pt), DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, tnew_dt); break; } case 2: { - umat_function_2(etot, Detot, sigma, Lt.slice(pt), DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, tnew_dt); + umat_function_finite(umat_name_py, etot, Detot, F0.slice(pt), F1.slice(pt), sigma, Lt.slice(pt), L.slice(pt), DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, tnew_dt); break; } - case 3: { - umat_function_3(etot, Detot, sigma, Lt.slice(pt), L.slice(pt), DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, tnew_dt); - break; - } - case 4: { - umat_function_4(umat_name_py, etot, Detot, F0.slice(pt), F1.slice(pt), sigma, Lt.slice(pt), L.slice(pt), DR.slice(pt), nprops, props, nstatev, statev, T, DT, Time, DTime, Wm(0), Wm(1), Wm(2), Wm(3), ndi, nshr, start, tnew_dt); - break; - } } } #ifdef _OPENMP diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 439b5125f..29d2a2ad1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -42,7 +42,6 @@ target_sources(simcoon PRIVATE Continuum_mechanics/Umat/Finite/generic_hyper_invariants.cpp Continuum_mechanics/Umat/Finite/generic_hyper_pstretch.cpp Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.cpp - Continuum_mechanics/Umat/Finite/mooney_rivlin.cpp Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp Continuum_mechanics/Umat/Finite/saint_venant.cpp @@ -72,7 +71,6 @@ target_sources(simcoon PRIVATE # Continuum Mechanics - UMAT Mechanical SMA Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp - Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.cpp Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp Continuum_mechanics/Umat/Mechanical/SMA/aniso_T.cpp diff --git a/src/Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp b/src/Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp index 8e2472e09..15f9ec397 100755 --- a/src/Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp +++ b/src/Continuum_mechanics/Umat/Finite/Neo_hookean_comp.cpp @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -45,9 +46,10 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_neo_hookean_comp(const vec &Etot, const vec &DEtot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) -{ +void umat_neo_hookean_comp(const string &umat_name, const vec &Etot, const vec &DEtot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +{ + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(statev); diff --git a/src/Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp b/src/Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp index ac703c234..040001e58 100755 --- a/src/Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp +++ b/src/Continuum_mechanics/Umat/Finite/Neo_hookean_incomp.cpp @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -43,9 +44,10 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_neo_hookean_incomp(const vec &etot, const vec &Detot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) -{ +void umat_neo_hookean_incomp(const string &umat_name, const vec &etot, const vec &Detot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +{ + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(statev); diff --git a/src/Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.cpp b/src/Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.cpp index c5f6148eb..f2826ed95 100755 --- a/src/Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.cpp +++ b/src/Continuum_mechanics/Umat/Finite/hypoelastic_orthotropic.cpp @@ -20,6 +20,7 @@ #include #include +#include #include #include #include @@ -38,9 +39,10 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_hypoelasticity_ortho(const vec &Etot, const vec &DEtot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) -{ +void umat_hypoelasticity_ortho(const string &umat_name, const vec &Etot, const vec &DEtot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +{ + UNUSED(umat_name); UNUSED(Etot); UNUSED(DR); UNUSED(nprops); @@ -94,12 +96,7 @@ void umat_hypoelasticity_ortho(const vec &Etot, const vec &DEtot, const mat &F0, vec DEel = DEtot - alpha*DT; sigma = el_pred(sigma_start, L, DEel); - if((solver_type == 0)||(solver_type==2)) { - Lt = L; - } - else if(solver_type == 1) { - sigma_in = zeros(6); - } + Lt = L; //Computation of the mechanical and thermal work quantities Wm += 0.5*sum((sigma_start+sigma)%DEtot); diff --git a/src/Continuum_mechanics/Umat/Finite/mooney_rivlin.cpp b/src/Continuum_mechanics/Umat/Finite/mooney_rivlin.cpp deleted file mode 100755 index f22455ae1..000000000 --- a/src/Continuum_mechanics/Umat/Finite/mooney_rivlin.cpp +++ /dev/null @@ -1,122 +0,0 @@ -/* This file is part of simcoon. - - simcoon is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - simcoon is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with simcoon. If not, see . - - */ - -///@file elastic_isotropic.cpp -///@brief User subroutine for Isotropic elastic materials in 3D case -///@version 1.0 - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace std; -using namespace arma; - -namespace simcoon{ - -///@brief The elastic UMAT requires 2 constants: -///@brief props[0] : Young modulus -///@brief props[1] : Poisson ratio -///@brief props[2] : CTE - -///@brief No statev is required for thermoelastic constitutive law - -void umat_mooney_rivlin(const vec &Etot, const vec &DEtot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) -{ - - UNUSED(nprops); - UNUSED(nstatev); - UNUSED(statev); - UNUSED(Time); - UNUSED(DTime); - UNUSED(nshr); - UNUSED(tnew_dt); - - double T_init = statev(0); - - //From the props to the material properties - double E = props(0); - double nu = props(1); - double alpha = props(2); - - double C_10 = E/(4.*(1+nu)); - double D_1 = 6.*(1-2.*nu)/E; - - ///@brief Initialization - if(start) - { - T_init = T; - sigma = zeros(6); - - Wm = 0.; - Wm_r = 0.; - Wm_ir = 0.; - Wm_d = 0.; - } - - vec sigma_start = sigma; - - //definition of the Right Cauchy-Green tensor - mat C = R_Cauchy_Green(F1); - - //Invariants of C - double I1 = trace(C); // pow(lambda_alpha(2),2.) + pow(lambda_alpha(1),2.) + pow(lambda_alpha(0),2.); //ascending order - double J = det(F1); //lambda(2)*lambda(1)*lambda(0) - double I1_bar = pow(J,-2./3.)*I1; - - double W = C_10*(I1_bar-3.) + (1./D_1)*pow(J-1.,2.); - - mat invC = inv(C); - mat I = eye(3,3); - - //Compute the PKII stress and then the Cauchy stress - mat S = (-2./3.)*C_10*I1_bar*invC + 2.*C_10*(1./J)*pow(J,-2./3.)*I + (2./D_1)*(J-1)*J*invC; - mat sigma_Cauchy = PKII2Cauchy(S, F1, J); - sigma = t2v_stress(sigma_Cauchy); - -// L = (-2./3.)*C_10*pow(J,-2./3.)*sym_dyadic(invC,I)+(2./9.)*C_10*I1_bar*sym_dyadic(invC,invC)-(2./3.)*C_10*I1_bar*dinvSdSsym(C) -// -(2./3.)*C_10*pow(J,-2./3.)*sym_dyadic(I,invC) -// +(1./D_1)*(J-1.)*J*sym_dyadic(invC,invC)+(2./D_1)*(J-1)*J*dinvSdSsym(C); - L = (-2./3.)*C_10*pow(J,-2./3.)*dyadic(invC,I)+(2./9.)*C_10*I1_bar*auto_dyadic(invC)-(2./3.)*C_10*I1_bar*dinvSdSsym(C) - -(2./3.)*C_10*pow(J,-2./3.)*dyadic(I,invC) - +(1./D_1)*(J-1.)*J*auto_dyadic(invC)+(2./D_1)*(J-1)*J*dinvSdSsym(C); - - - if((solver_type == 0)||(solver_type==2)) { - Lt = L; - } - - //Computation of the mechanical and thermal work quantities - /* - Wm += 0.5*sum((sigma_start+sigma)%DEtot); - Wm_r += 0.5*sum((sigma_start+sigma)%DEtot); - Wm_ir += 0.; - Wm_d += 0.; - */ - - statev(0) = T_init; -} - -} //namespace simcoon diff --git a/src/Continuum_mechanics/Umat/Finite/saint_venant.cpp b/src/Continuum_mechanics/Umat/Finite/saint_venant.cpp index ca2aaa0fe..bb6c5ea27 100755 --- a/src/Continuum_mechanics/Umat/Finite/saint_venant.cpp +++ b/src/Continuum_mechanics/Umat/Finite/saint_venant.cpp @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -45,9 +46,10 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_saint_venant(const vec &etot, const vec &Detot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) -{ +void umat_saint_venant(const string &umat_name, const vec &etot, const vec &Detot, const mat &F0, const mat &F1, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +{ + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(statev); diff --git a/src/Continuum_mechanics/Umat/Mechanical/Combined/Prony_Nfast_Plastic.cpp b/src/Continuum_mechanics/Umat/Mechanical/Combined/Prony_Nfast_Plastic.cpp index 68dbca7a3..407492216 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Combined/Prony_Nfast_Plastic.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Combined/Prony_Nfast_Plastic.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -30,9 +31,10 @@ using namespace arma; namespace simcoon { -void umat_prony_Nfast_plastic(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +void umat_prony_Nfast_plastic(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -67,6 +69,7 @@ void umat_prony_Nfast_plastic(const vec &Etot, const vec &DEtot, vec &sigma, mat //Define the viscoelastic stiffness mat L0 = L_iso(E0, nu0, "Enu"); + L = L0; mat M0 = M_iso(E0, nu0, "Enu"); ///@brief Temperature initialization double T_init = statev(0); diff --git a/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_LLD_0.cpp b/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_LLD_0.cpp index 688c82f4c..0c4822169 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_LLD_0.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_LLD_0.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -46,8 +47,9 @@ namespace simcoon { ///@param props(5) : lambdaD Damage evolution parameter lambda ///@param props(6) : deltaD Damage evolution parameter delta -void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) { - +void umat_damage_LLD_0(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -480,8 +482,7 @@ void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, m Eel = Etot + DEtot - alpha*(T+DT-Tinit) - EP; sigma = el_pred(L_tilde, Eel, ndi); - if((solver_type == 0)||(solver_type == 2)) { - + { //Tangent modulus mat B = L*inv(L_tilde); //stress "localization factor" in damage @@ -619,9 +620,6 @@ void umat_damage_LLD_0(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, m Lt = L_tilde - (kappa_j[0]*P_epsilon[0].t() + kappa_j[1]*P_epsilon[1].t() + kappa_j[2]*P_epsilon[2].t()); } - else if(solver_type == 1) { - sigma_in = -L*(Etot - Eel); - } statev(0) = Tinit; statev(1) = d_22; diff --git a/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.cpp b/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.cpp index bb6bac4e4..b259c6ff4 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Damage/damage_weibull.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include @@ -40,8 +41,9 @@ namespace simcoon { ///@param props(5) : lambdaD Damage evolution parameter lambda ///@param props(6) : deltaD Damage evolution parameter delta -void umat_damage_weibull(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) { - +void umat_damage_weibull(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -49,8 +51,6 @@ void umat_damage_weibull(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, UNUSED(nshr); UNUSED(tnew_dt); UNUSED(DR); - UNUSED(sigma_in); - UNUSED(solver_type); //From the props to the material properties double E = props(0); diff --git a/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_isotropic.cpp b/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_isotropic.cpp index 6bd15bd7c..46cb295bd 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_isotropic.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_isotropic.cpp @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -38,9 +39,10 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_elasticity_iso(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) -{ +void umat_elasticity_iso(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +{ + UNUSED(umat_name); UNUSED(Etot); UNUSED(DR); UNUSED(nprops); @@ -50,9 +52,9 @@ void umat_elasticity_iso(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, UNUSED(DTime); UNUSED(nshr); UNUSED(tnew_dt); - + double T_init = statev(0); - + //From the props to the material properties double E = props(0); double nu = props(1); @@ -60,38 +62,33 @@ void umat_elasticity_iso(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, //Elastic stiffness tensor L = L_iso(E, nu, "Enu"); - + ///@brief Initialization if(start) { T_init = T; sigma = zeros(6); - + Wm = 0.; Wm_r = 0.; Wm_ir = 0.; Wm_d = 0.; } - + vec sigma_start = sigma; - - //Compute the elastic strain and the related stress + + //Compute the elastic strain and the related stress vec Eel = Etot + DEtot - alpha*(T+DT-T_init); sigma = el_pred(L, Eel, ndi); - - if((solver_type == 0)||(solver_type==2)) { - Lt = L; - } - else if(solver_type == 1) { - sigma_in = zeros(6); - } - + + Lt = L; + //Computation of the mechanical and thermal work quantities Wm += 0.5*sum((sigma_start+sigma)%DEtot); Wm_r += 0.5*sum((sigma_start+sigma)%DEtot); Wm_ir += 0.; Wm_d += 0.; - + statev(0) = T_init; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_orthotropic.cpp b/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_orthotropic.cpp index aad5ac1a4..e89c19152 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_orthotropic.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_orthotropic.cpp @@ -20,6 +20,7 @@ #include #include +#include #include #include #include @@ -38,9 +39,10 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_elasticity_ortho(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) -{ +void umat_elasticity_ortho(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +{ + UNUSED(umat_name); UNUSED(Etot); UNUSED(DR); UNUSED(nprops); @@ -50,9 +52,9 @@ void umat_elasticity_ortho(const vec &Etot, const vec &DEtot, vec &sigma, mat &L UNUSED(DTime); UNUSED(nshr); UNUSED(tnew_dt); - + double T_init = statev(0); - + //From the props to the material properties double Ex = props(0); double Ey = props(1); @@ -68,45 +70,40 @@ void umat_elasticity_ortho(const vec &Etot, const vec &DEtot, vec &sigma, mat &L double alphaz = props(11); //Elastic stiffness tensor - L = L_ortho(Ex,Ey,Ez,nuxy,nuxz,nuyz,Gxy,Gxz,Gyz, "EnuG"); - + L = L_ortho(Ex,Ey,Ez,nuxy,nuxz,nuyz,Gxy,Gxz,Gyz, "EnuG"); + ///@brief Initialization if(start) { T_init = T; sigma = zeros(6); - + Wm = 0.; Wm_r = 0.; Wm_ir = 0.; Wm_d = 0.; } - + vec sigma_start = sigma; - + //definition of the CTE tensor vec alpha = zeros(6); alpha(0) = alphax; alpha(1) = alphay; alpha(2) = alphaz; - - //Compute the elastic strain and the related stress + + //Compute the elastic strain and the related stress vec Eel = Etot + DEtot - alpha*(T+DT-T_init); sigma = el_pred(L, Eel, ndi); - - if((solver_type == 0)||(solver_type==2)) { - Lt = L; - } - else if(solver_type == 1) { - sigma_in = zeros(6); - } - + + Lt = L; + //Computation of the mechanical and thermal work quantities Wm += 0.5*sum((sigma_start+sigma)%DEtot); Wm_r += 0.5*sum((sigma_start+sigma)%DEtot); Wm_ir += 0.; Wm_d += 0.; - + statev(0) = T_init; } diff --git a/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_transverse_isotropic.cpp b/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_transverse_isotropic.cpp index 18c0eed73..2c9baa7ea 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_transverse_isotropic.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Elasticity/elastic_transverse_isotropic.cpp @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -43,9 +44,10 @@ namespace simcoon{ ///@brief No statev is required for thermoelastic constitutive law -void umat_elasticity_trans_iso(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) -{ +void umat_elasticity_trans_iso(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +{ + UNUSED(umat_name); UNUSED(Etot); UNUSED(DR); UNUSED(nprops); @@ -55,9 +57,9 @@ void umat_elasticity_trans_iso(const vec &Etot, const vec &DEtot, vec &sigma, ma UNUSED(DTime); UNUSED(nshr); UNUSED(tnew_dt); - + double T_init = statev(0); - + //From the props to the material properties double axis = props(0); double EL = props(1); @@ -69,47 +71,42 @@ void umat_elasticity_trans_iso(const vec &Etot, const vec &DEtot, vec &sigma, ma double alphaT = props(7); //Elastic stiffness tensor - L = L_isotrans(EL, ET, nuTL, nuTT, GLT, axis); - + L = L_isotrans(EL, ET, nuTL, nuTT, GLT, axis); + ///@brief Initialization if(start) { T_init = T; sigma = zeros(6); - + Wm = 0.; Wm_r = 0.; Wm_ir = 0.; Wm_d = 0.; } - + vec sigma_start = sigma; //definition of the CTE tensor vec alpha = zeros(6); alpha = alphaT*Ith(); alpha(axis-1) += alphaL-alphaT; - + ///Elastic prediction - Accounting for the thermal prediction - //Compute the elastic strain and the related stress + //Compute the elastic strain and the related stress vec Eel = Etot + DEtot - alpha*(T+DT-T_init); sigma = el_pred(L, Eel, ndi); - - if((solver_type == 0)||(solver_type==2)) { - Lt = L; - } - else if(solver_type == 1) { - sigma_in = zeros(6); - } - + + Lt = L; + //Computation of the mechanical and thermal work quantities Wm += 0.5*sum((sigma_start+sigma)%DEtot); Wm_r += 0.5*sum((sigma_start+sigma)%DEtot); Wm_ir += 0.; Wm_d += 0.; - - statev(0) = T_init; + + statev(0) = T_init; } } //namespace simcoon diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.cpp index c9ba568a9..cf2c18be3 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -64,9 +65,10 @@ namespace simcoon { ///@brief statev[13] : Backstress 11: X(1,2) -void umat_ani_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) +void umat_ani_chaboche_CCP(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -287,46 +289,40 @@ void umat_ani_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &L vec Da_1 = a_1 - a_1start; vec Da_2 = a_2 - a_2start; - if((solver_type == 0)||(solver_type==2)) { - - //Computation of the tangent modulus - mat Bhat = zeros(1, 1); - Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); - - vec op = zeros(1); - mat delta = eye(1,1); - - for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) - op(i) = 1.; - } - - mat Bbar = zeros(1,1); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); - } - } - - mat invBbar = zeros(1, 1); - mat invBhat = zeros(1, 1); - invBbar = inv(Bbar); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - invBhat(i, j) = op(i)*op(j)*invBbar(i, j); - } - } - - std::vector P_epsilon(1); - P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); - std::vector P_theta(1); - P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); - - Lt = L - (kappa_j[0]*P_epsilon[0].t()); - } - else if(solver_type == 1) { - sigma_in = -L*EP; + //Computation of the tangent modulus + mat Bhat = zeros(1, 1); + Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); + + vec op = zeros(1); + mat delta = eye(1,1); + + for (int i=0; i<1; i++) { + if(Ds_j[i] > sim_iota) + op(i) = 1.; } + + mat Bbar = zeros(1,1); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); + } + } + + mat invBbar = zeros(1, 1); + mat invBhat = zeros(1, 1); + invBbar = inv(Bbar); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + invBhat(i, j) = op(i)*op(j)*invBbar(i, j); + } + } + + std::vector P_epsilon(1); + P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); + std::vector P_theta(1); + P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); + + Lt = L - (kappa_j[0]*P_epsilon[0].t()); double A_p = -Hp; vec A_a1 = -X_1; diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.cpp index 28452ebd0..2d435e63f 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -64,9 +65,10 @@ namespace simcoon { ///@brief statev[13] : Backstress 11: X(1,2) -void umat_dfa_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) +void umat_dfa_chaboche_CCP(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -285,46 +287,40 @@ void umat_dfa_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &L vec Da_1 = a_1 - a_1start; vec Da_2 = a_2 - a_2start; - if((solver_type == 0)||(solver_type==2)) { - - //Computation of the tangent modulus - mat Bhat = zeros(1, 1); - Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); - - vec op = zeros(1); - mat delta = eye(1,1); - - for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) - op(i) = 1.; - } - - mat Bbar = zeros(1,1); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); - } - } - - mat invBbar = zeros(1, 1); - mat invBhat = zeros(1, 1); - invBbar = inv(Bbar); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - invBhat(i, j) = op(i)*op(j)*invBbar(i, j); - } - } - - std::vector P_epsilon(1); - P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); - std::vector P_theta(1); - P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); - - Lt = L - (kappa_j[0]*P_epsilon[0].t()); - } - else if(solver_type == 1) { - sigma_in = -L*EP; + //Computation of the tangent modulus + mat Bhat = zeros(1, 1); + Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); + + vec op = zeros(1); + mat delta = eye(1,1); + + for (int i=0; i<1; i++) { + if(Ds_j[i] > sim_iota) + op(i) = 1.; } + + mat Bbar = zeros(1,1); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); + } + } + + mat invBbar = zeros(1, 1); + mat invBhat = zeros(1, 1); + invBbar = inv(Bbar); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + invBhat(i, j) = op(i)*op(j)*invBbar(i, j); + } + } + + std::vector P_epsilon(1); + P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); + std::vector P_theta(1); + P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); + + Lt = L - (kappa_j[0]*P_epsilon[0].t()); double A_p = -Hp; vec A_a1 = -X_1; diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.cpp index d5a5b32f5..95952823b 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Generic_chaboche_ccp.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -64,9 +65,10 @@ namespace simcoon { ///@brief statev[13] : Backstress 11: X(1,2) -void umat_generic_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) +void umat_generic_chaboche_CCP(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -434,46 +436,40 @@ void umat_generic_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, ma Da_N[i] = a_N[i] - a_N_start[i]; } - if((solver_type == 0)||(solver_type==2)) { - - //Computation of the tangent modulus - mat Bhat = zeros(1, 1); - Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); - - vec op = zeros(1); - mat delta = eye(1,1); - - for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) - op(i) = 1.; - } - - mat Bbar = zeros(1,1); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); - } - } - - mat invBbar = zeros(1, 1); - mat invBhat = zeros(1, 1); - invBbar = inv(Bbar); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - invBhat(i, j) = op(i)*op(j)*invBbar(i, j); - } - } - - std::vector P_epsilon(1); - P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); - std::vector P_theta(1); - P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); - - Lt = L - (kappa_j[0]*P_epsilon[0].t()); - } - else if(solver_type == 1) { - sigma_in = -L*EP; + //Computation of the tangent modulus + mat Bhat = zeros(1, 1); + Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); + + vec op = zeros(1); + mat delta = eye(1,1); + + for (int i=0; i<1; i++) { + if(Ds_j[i] > sim_iota) + op(i) = 1.; } + + mat Bbar = zeros(1,1); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); + } + } + + mat invBbar = zeros(1, 1); + mat invBhat = zeros(1, 1); + invBbar = inv(Bbar); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + invBhat(i, j) = op(i)*op(j)*invBbar(i, j); + } + } + + std::vector P_epsilon(1); + P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); + std::vector P_theta(1); + P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); + + Lt = L - (kappa_j[0]*P_epsilon[0].t()); double A_p = -Hp; diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.cpp index e15d0c295..d9873ab8f 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_chaboche_ccp.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -64,9 +65,10 @@ namespace simcoon { ///@brief statev[13] : Backstress 11: X(1,2) -void umat_hill_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) +void umat_hill_chaboche_CCP(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -284,46 +286,40 @@ void umat_hill_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat & vec Da_1 = a_1 - a_1start; vec Da_2 = a_2 - a_2start; - if((solver_type == 0)||(solver_type==2)) { - - //Computation of the tangent modulus - mat Bhat = zeros(1, 1); - Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); - - vec op = zeros(1); - mat delta = eye(1,1); - - for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) - op(i) = 1.; - } - - mat Bbar = zeros(1,1); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); - } - } - - mat invBbar = zeros(1, 1); - mat invBhat = zeros(1, 1); - invBbar = inv(Bbar); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - invBhat(i, j) = op(i)*op(j)*invBbar(i, j); - } - } - - std::vector P_epsilon(1); - P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); - std::vector P_theta(1); - P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); - - Lt = L - (kappa_j[0]*P_epsilon[0].t()); - } - else if(solver_type == 1) { - sigma_in = -L*EP; + //Computation of the tangent modulus + mat Bhat = zeros(1, 1); + Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); + + vec op = zeros(1); + mat delta = eye(1,1); + + for (int i=0; i<1; i++) { + if(Ds_j[i] > sim_iota) + op(i) = 1.; } + + mat Bbar = zeros(1,1); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); + } + } + + mat invBbar = zeros(1, 1); + mat invBhat = zeros(1, 1); + invBbar = inv(Bbar); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + invBhat(i, j) = op(i)*op(j)*invBbar(i, j); + } + } + + std::vector P_epsilon(1); + P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); + std::vector P_theta(1); + P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); + + Lt = L - (kappa_j[0]*P_epsilon[0].t()); double A_p = -Hp; vec A_a1 = -X_1; diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.cpp index 3865fac2f..175556358 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -80,9 +81,10 @@ namespace simcoon { A_theta = */ -void umat_plasticity_hill_isoh_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +void umat_plasticity_hill_isoh_CCP(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -111,7 +113,7 @@ void umat_plasticity_hill_isoh_CCP(const vec &Etot, const vec &DEtot, vec &sigma vec alpha = alpha_iso*Ith(); //Elastic stiffness tensors - mat L = L_iso(E, nu, "Enu"); + L = L_iso(E, nu, "Enu"); mat M = M_iso(E, nu, "Enu"); ///@brief Temperature initialization diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.cpp index 998de2678..de2b6c138 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/Hill_isoh_Nfast.cpp @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -47,9 +48,10 @@ using namespace arma; namespace simcoon { -void umat_plasticity_hill_isoh_CCP_N(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +void umat_plasticity_hill_isoh_CCP_N(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -94,7 +96,7 @@ void umat_plasticity_hill_isoh_CCP_N(const vec &Etot, const vec &DEtot, vec &sig vec alpha = alpha_iso*Ith(); //Define the elastic stiffness - mat L = L_iso(E, nu, "Enu"); + L = L_iso(E, nu, "Enu"); mat M = M_iso(E, nu, "Enu"); ///@brief Temperature initialization double T_init = statev(0); diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.cpp index 2f61984a6..cc160a0e3 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_chaboche_ccp.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -63,9 +64,10 @@ namespace simcoon { ///@brief statev[13] : Backstress 11: X(1,2) -void umat_plasticity_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) +void umat_plasticity_chaboche_CCP(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -272,46 +274,40 @@ void umat_plasticity_chaboche_CCP(const vec &Etot, const vec &DEtot, vec &sigma, vec Da_1 = a_1 - a_1start; vec Da_2 = a_2 - a_2start; - if((solver_type == 0)||(solver_type==2)) { - - //Computation of the tangent modulus - mat Bhat = zeros(1, 1); - Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); - - vec op = zeros(1); - mat delta = eye(1,1); - - for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) - op(i) = 1.; - } - - mat Bbar = zeros(1,1); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); - } - } - - mat invBbar = zeros(1, 1); - mat invBhat = zeros(1, 1); - invBbar = inv(Bbar); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - invBhat(i, j) = op(i)*op(j)*invBbar(i, j); - } - } - - std::vector P_epsilon(1); - P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); - std::vector P_theta(1); - P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); - - Lt = L - (kappa_j[0]*P_epsilon[0].t()); - } - else if(solver_type == 1) { - sigma_in = -L*EP; + //Computation of the tangent modulus + mat Bhat = zeros(1, 1); + Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); + + vec op = zeros(1); + mat delta = eye(1,1); + + for (int i=0; i<1; i++) { + if(Ds_j[i] > sim_iota) + op(i) = 1.; } + + mat Bbar = zeros(1,1); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); + } + } + + mat invBbar = zeros(1, 1); + mat invBhat = zeros(1, 1); + invBbar = inv(Bbar); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + invBhat(i, j) = op(i)*op(j)*invBbar(i, j); + } + } + + std::vector P_epsilon(1); + P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); + std::vector P_theta(1); + P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); + + Lt = L - (kappa_j[0]*P_epsilon[0].t()); double A_p = -Hp; vec A_a1 = -X_1; diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.cpp index fa4238780..cc02e983c 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_isotropic_ccp.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -81,9 +82,10 @@ namespace simcoon { A_theta = */ -void umat_plasticity_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) +void umat_plasticity_iso_CCP(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -221,47 +223,41 @@ void umat_plasticity_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat vec DEP = EP - EP_start; double Dp = Ds_j[0]; - if((solver_type == 0)||(solver_type==2)) { - - //Computation of the tangent modulus - mat Bhat = zeros(1, 1); - Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); - - vec op = zeros(1); - mat delta = eye(1,1); - - for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) - op(i) = 1.; - } - - mat Bbar = zeros(1,1); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); - } - } - - mat invBbar = zeros(1, 1); - mat invBhat = zeros(1, 1); - invBbar = inv(Bbar); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - invBhat(i, j) = op(i)*op(j)*invBbar(i, j); - } + //Computation of the tangent modulus + mat Bhat = zeros(1, 1); + Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); + + vec op = zeros(1); + mat delta = eye(1,1); + + for (int i=0; i<1; i++) { + if(Ds_j[i] > sim_iota) + op(i) = 1.; + } + + mat Bbar = zeros(1,1); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); } - - std::vector P_epsilon(1); - P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); - std::vector P_theta(1); - P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); - - Lt = L - (kappa_j[0]*P_epsilon[0].t()); } - else if(solver_type == 1) { - sigma_in = -L*EP; + + mat invBbar = zeros(1, 1); + mat invBhat = zeros(1, 1); + invBbar = inv(Bbar); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + invBhat(i, j) = op(i)*op(j)*invBbar(i, j); + } } + std::vector P_epsilon(1); + P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); + std::vector P_theta(1); + P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); + + Lt = L - (kappa_j[0]*P_epsilon[0].t()); + double A_p = -Hp; double Dgamma_loc = 0.5*sum((sigma_start+sigma)%DEP) + 0.5*(A_p_start + A_p)*Dp; diff --git a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.cpp b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.cpp index 713660794..a8e8345b6 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Plasticity/plastic_kin_iso_ccp.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -63,9 +64,10 @@ namespace simcoon { ///@brief statev[13] : Backstress 11: X(1,2) -void umat_plasticity_kin_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) +void umat_plasticity_kin_iso_CCP(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -229,46 +231,40 @@ void umat_plasticity_kin_iso_CCP(const vec &Etot, const vec &DEtot, vec &sigma, double Dp = Ds_j[0]; vec Da = a - a_start; - if((solver_type == 0)||(solver_type==2)) { - - //Computation of the tangent modulus - mat Bhat = zeros(1, 1); - Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); - - vec op = zeros(1); - mat delta = eye(1,1); - - for (int i=0; i<1; i++) { - if(Ds_j[i] > sim_iota) - op(i) = 1.; - } - - mat Bbar = zeros(1,1); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); - } - } - - mat invBbar = zeros(1, 1); - mat invBhat = zeros(1, 1); - invBbar = inv(Bbar); - for (int i = 0; i < 1; i++) { - for (int j = 0; j < 1; j++) { - invBhat(i, j) = op(i)*op(j)*invBbar(i, j); - } - } - - std::vector P_epsilon(1); - P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); - std::vector P_theta(1); - P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); - - Lt = L - (kappa_j[0]*P_epsilon[0].t()); - } - else if(solver_type == 1) { - sigma_in = -L*EP; + //Computation of the tangent modulus + mat Bhat = zeros(1, 1); + Bhat(0, 0) = sum(dPhidsigma%kappa_j[0]) - K(0,0); + + vec op = zeros(1); + mat delta = eye(1,1); + + for (int i=0; i<1; i++) { + if(Ds_j[i] > sim_iota) + op(i) = 1.; } + + mat Bbar = zeros(1,1); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + Bbar(i, j) = op(i)*op(j)*Bhat(i, j) + delta(i,j)*(1-op(i)*op(j)); + } + } + + mat invBbar = zeros(1, 1); + mat invBhat = zeros(1, 1); + invBbar = inv(Bbar); + for (int i = 0; i < 1; i++) { + for (int j = 0; j < 1; j++) { + invBhat(i, j) = op(i)*op(j)*invBbar(i, j); + } + } + + std::vector P_epsilon(1); + P_epsilon[0] = invBhat(0, 0)*(L*dPhidsigma); + std::vector P_theta(1); + P_theta[0] = dPhidtheta - sum(dPhidsigma%(L*alpha)); + + Lt = L - (kappa_j[0]*P_epsilon[0].t()); double A_p = -Hp; vec A_a = -X; diff --git a/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp b/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp index 4af00c628..e5c4a6eb1 100644 --- a/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono.cpp @@ -16,11 +16,12 @@ */ ///@file sma_mono.cpp -///@brief User subroutine for monocrystalline SMA +///@brief User subroutine for monocrystalline SMA with multiple elastic symmetry options ///@author Chemisky #include #include +#include #include #include @@ -35,23 +36,20 @@ using namespace arma; namespace simcoon { -///@brief The elastic UMAT requires 17 constants: -///@brief props(0) : E -///@brief props(1) : nu -///@brief props(2) : alpha_iso -///@brief props(3) : b It is the slope that correspond to g -///@brief props(4) : g the shear strain -///@brief props(5) : Ms Martensitic start temperature -///@brief props(6) : Af Austenitic finish temperature -///@brief props(7) : nvariants The number of martensitic variants -///@brief props(8) : c_lambda0 -///@brief props(9) : p0_lambda0 -///@brief props(10) : n_lambda0 -///@brief props(11) : alpha_lambda0 -///@brief props(12) : c_lambda1 -///@brief props(13) : p0_lambda1 -///@brief props(14) : n_lambda1 -///@brief props(15) : alpha_lambda1 +///@brief Unified SMA monocrystal UMAT supporting multiple elastic symmetries +///@brief The elastic symmetry is determined by umat_name: +///@brief SMAMO (isotropic): props(0)=E, props(1)=nu, then common params from props(2) +///@brief SMAMC (cubic): props(0)=C11, props(1)=C12, props(2)=C44, then common params from props(3) +///@brief +///@brief Common parameters (offset depends on symmetry): +///@brief alpha_iso : Coefficient of thermal expansion +///@brief b : Slope parameter (corresponds to g) +///@brief g : Shear strain magnitude +///@brief Ms : Martensite start temperature +///@brief Af : Austenite finish temperature +///@brief nvariants : Number of martensite variants +///@brief c_lambda0, p0_lambda0, n_lambda0, alpha_lambda0 : Lagrange parameters (variant) +///@brief c_lambda1, p0_lambda1, n_lambda1, alpha_lambda1 : Lagrange parameters (total) /*void Amortissement(double &xi, const double &dxi, const double &xi_start) { @@ -71,40 +69,77 @@ namespace simcoon { } }*/ -void umat_sma_mono(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - +void umat_sma_mono(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { + UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); UNUSED(DTime); UNUSED(nshr); UNUSED(tnew_dt); - + // ###################### Props ################################# //From the props to the material properties - double E = props(0); - double nu = props(1); - double alpha_iso = props(2); - double b = props(3); - double g = props(4); - double Ms = props(5); - double Af = props(6); - int nvariants = props(7); - double c_lambda0 = props(8); - double p0_lambda0 = props(9); - double n_lambda0 = props(10); - double alpha_lambda0 = props(11); - double c_lambda1 = props(12); - double p0_lambda1 = props(13); - double n_lambda1 = props(14); - double alpha_lambda1 = props(15); - + // Property offset depends on elastic symmetry type + int offset = 0; + double alpha_iso = 0.; + double b = 0.; + double g = 0.; + double Ms = 0.; + double Af = 0.; + int nvariants = 0; + double c_lambda0 = 0.; + double p0_lambda0 = 0.; + double n_lambda0 = 0.; + double alpha_lambda0 = 0.; + double c_lambda1 = 0.; + double p0_lambda1 = 0.; + double n_lambda1 = 0.; + double alpha_lambda1 = 0.; + + // Build elastic stiffness tensor based on umat_name + if (umat_name == "SMAMO") { + // Isotropic elasticity: E, nu + double E = props(0); + double nu = props(1); + L = L_iso(E, nu, "Enu"); + offset = 2; + } + else if (umat_name == "SMAMC") { + // Cubic elasticity: C11, C12, C44 + double C11 = props(0); + double C12 = props(1); + double C44 = props(2); + L = L_cubic(C11, C12, C44, "Cii"); + offset = 3; + } + else { + cout << "Error: Unknown umat_name in umat_sma_mono: " << umat_name << "\n"; + exit(0); + } + + // Extract common parameters using offset + alpha_iso = props(offset); + b = props(offset + 1); + g = props(offset + 2); + Ms = props(offset + 3); + Af = props(offset + 4); + nvariants = int(props(offset + 5)); + c_lambda0 = props(offset + 6); + p0_lambda0 = props(offset + 7); + n_lambda0 = props(offset + 8); + alpha_lambda0 = props(offset + 9); + c_lambda1 = props(offset + 10); + p0_lambda1 = props(offset + 11); + n_lambda1 = props(offset + 12); + alpha_lambda1 = props(offset + 13); + double roS0 = -1.*b*g; double romu0 = 0.5*roS0*(Ms + Af); double Y0 = 0.5*roS0*(Ms - Af); - + //definition of the CTE tensor - vec alpha = alpha_iso*Ith(); + vec alpha = alpha_iso*Ith(); std::string data_path= std::getenv("SIMCOON_DATA_PATH") ? std::getenv("SIMCOON_DATA_PATH") : "data" ; mat Hnm = zeros(nvariants, nvariants); Hnm.load(data_path+"/Hnm.inp", raw_ascii); @@ -157,9 +192,8 @@ void umat_sma_mono(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat & //Rotation of internal variables (tensors) ET = rotate_strain(ET, DR); - //Elstic stiffness tensor - L = L_iso(E, nu, "Enu"); - + // Note: Elastic stiffness tensor L is already constructed based on umat_name + ///@brief Initialization if(start) { diff --git a/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.cpp b/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.cpp deleted file mode 100644 index 541b5c4a7..000000000 --- a/src/Continuum_mechanics/Umat/Mechanical/SMA/SMA_mono_cubic.cpp +++ /dev/null @@ -1,445 +0,0 @@ -/* This file is part of simcoon. - - simcoon is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - simcoon is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with simcoon. If not, see . - - */ - -///@file sma_mono.cpp -///@brief User subroutine for monocrystalline SMA -///@author Chemisky - -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -using namespace std; -using namespace arma; - -namespace simcoon { - -///@brief The elastic UMAT requires 18 constants: -///@brief props(0) : C11 -///@brief props(1) : C12 -///@brief props(2) : C44 -///@brief props(3) : alpha_iso -///@brief props(4) : b It is the slope that correspond to g -///@brief props(5) : g the shear strain -///@brief props(6) : Ms Martensitic start temperature -///@brief props(7) : Af Austenitic finish temperature -///@brief props(8) : nvariants The number of martensitic variants -///@brief props(9) : c_lambda0 -///@brief props(10) : p0_lambda0 -///@brief props(11) : n_lambda0 -///@brief props(12) : alpha_lambda0 -///@brief props(13) : c_lambda1 -///@brief props(14) : p0_lambda1 -///@brief props(15) : n_lambda1 -///@brief props(16) : alpha_lambda1 - -///@brief No statev is required for thermoelastic constitutive law - -/*void Amortissement(double &xi, const double &dxi, const double &xi_start) -{ - - double limit1 = 1. - limit; - double limit2 = limit; - - - if (xi >= limit1) - { - xi = (xi - dxi + limit1)*0.5; - } - - if (xi <= limit2) - { - xi = (xi - dxi + limit2)*0.5; - } -}*/ - -void umat_sma_mono_cubic(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - - UNUSED(nprops); - UNUSED(nstatev); - UNUSED(Time); - UNUSED(DTime); - UNUSED(nshr); - UNUSED(tnew_dt); - - // ###################### Props ################################# - //From the props to the material properties - double C11 = props(0); - double C12 = props(1); - double C44 = props(2); - double alpha_iso = props(3); - double b = props(4); - double g = props(5); - double Ms = props(6); - double Af = props(7); - int nvariants = props(8); - double c_lambda0 = props(9); - double p0_lambda0 = props(10); - double n_lambda0 = props(11); - double alpha_lambda0 = props(12); - double c_lambda1 = props(13); - double p0_lambda1 = props(14); - double n_lambda1 = props(15); - double alpha_lambda1 = props(16); - - double roS0 = -1.*b*g; - double romu0 = 0.5*roS0*(Ms + Af); - double Y0 = 0.5*roS0*(Ms - Af); - - //definition of the CTE tensor - vec alpha = alpha_iso*Ith(); - - std::string data_path= std::getenv("SIMCOON_DATA_PATH") ? std::getenv("SIMCOON_DATA_PATH") : "data" ; - mat Hnm = zeros(nvariants, nvariants); - Hnm.load(data_path+"/Hnm.inp", raw_ascii); - - // ###################### Statev ################################# - - ///@brief Temperature initialization - double T_init = statev(0); - - std::vector var(nvariants); - - ///@brief Properties of the variants, use "test.dat" to specify the parameters (for now) - ifstream paramvariant; - paramvariant.open(data_path+"/variant.inp", ios::in); - if(paramvariant) { - string chaine1; - for(int i=0; i> chaine1 >> var[i].n(0) >> var[i].n(1) >> var[i].n(2) >> var[i].m(0) >> var[i].m(1) >> var[i].m(2); - var[i].build(g); - } - } - else { - cout << "Error: cannot open .dat file \n"; - } - paramvariant.close(); - - vec xin(nvariants); - vec xin_start(nvariants); - - for(int i=0; i transfo_actif(24); - - ///Elastic prediction - Accounting for the thermal prediction - vec Eel = Etot + DEtot - alpha*(T+DT-T_init) - ET; - sigma = el_pred(L, Eel, ndi); - - //Need to define the thermodynamic parameters: - vec lambda0 = zeros(nvariants); - - double lambda1 = 0.; - - mat dlambda0dxin = zeros(nvariants, nvariants); - - double dlambda1dxin = 0.; - - vec Hnmxim = Hnm*xin; - vec Hnmxim_start = Hnm*xin_start; - - lambda1 = lagrange_pow_1(xi, c_lambda1, p0_lambda1, n_lambda1, alpha_lambda1); - - //Set the Lagrange multiplliers due to Physical limitations for each Phi - for(int i=0; i 0)&&(PhiF(i) > PhiF_start(i))) { - transfo_actif[i] = 1; - } - else if((PhiR(i) < 0)&&(PhiR(i) < PhiR_start(i))) { - transfo_actif[i] = -1; - }*/ - - if(PhiF(i) > 0) { - transfo_actif[i] = 1; - } - else if(PhiR(i) < 0) { - transfo_actif[i] = -1; - } - - } - - //Collect the number of potential transforming systems - //Note : It is possible that some systems have a reverse transformation while others have a forward one .. - int nactive = 0; - for(int i=0; i 0) { - nactive++; - } - } - - vec dxin = zeros(nactive); - vec Phi = zeros(nactive); - mat dPhidxi = zeros(nactive, nactive); - mat dPhidsigma = zeros(nactive, 6); - mat lambda = zeros(6, nactive); - - std::vector active; - std::vector tactive; - - double sumPhi = 0.; - for(int i=0; i 0) { - sumPhi += fabs(PhiF(i)); - active.push_back(i); - tactive.push_back(1); - } - else if (transfo_actif[i] < 0) { - sumPhi += fabs(PhiR(i)); - active.push_back(i); - tactive.push_back(-1); - } - } - - int compteur = 0; - - if(nactive > 0) - { - - for(compteur = 0; ((compteur < maxiter_umat) && (sumPhi/Y0 > precision_umat)); compteur++) - { - dPhidxi = zeros(nactive, nactive); - - lambda1 = lagrange_pow_1(xi, c_lambda1, p0_lambda1, n_lambda1, alpha_lambda1); - dlambda1dxin = dlagrange_pow_1(xi, c_lambda1, p0_lambda1, n_lambda1, alpha_lambda1); - - for(int i=0; i 0.) - dxin = -1.*inv(dPhidxi)*Phi; - else { - //pnewdt = 0.1; -// Phi = zeros(nactive); - dxin = zeros(nactive); - } - - dxi = 0.; - for(int i=0; i 1.E-8) { - - xi_path = xi - xi_start; - - for(int i=0; i 0.)&&(fabs(xi - xi_start) > 1.E-5)) { - if(fabs(det(dPhidxi)) > 0.) { - Lt = L + ((L*lambda)*(inv(dPhidxi))*(dPhidsigma*L)); - } - else - Lt = L; - - vec eigval = eig_sym(Lt); - - ///Note : To use with Self-consistent micromechanics only! -/* while(eigval(0) < 0.) { - Lt = 0.99*Lt + 0.01*Lt_eff; - eigval = eig_sym(Lt); - } */ - - } - else { - - Eel = Etot + DEtot - alpha*(T+DT-T_init) - ET; - sigma = el_pred(L, Eel, ndi); - - ///Computation of the tangent modulus - Lt = L; - } - -/* if (compteur == maxitNewton_mono) { - tnew_dt*=0.2; - } - else if (compteur <= incNewton_mono ) { - tnew_dt*= 2.; - }*/ - - statev(0) = T_init; - - for(int i=0; i #include #include +#include #include #include #include @@ -39,8 +40,9 @@ using namespace arma; namespace simcoon{ -void umat_sma_aniso_T(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - +void umat_sma_aniso_T(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -196,7 +198,7 @@ void umat_sma_aniso_T(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, co mat M_A = M_iso(K_A, mu_A, "Kmu"); mat M_M = M_iso(K_M, mu_M, "Kmu"); mat M = M_iso(K_eff, mu_eff, "Kmu"); - mat L = L_iso(K_eff, mu_eff, "Kmu"); + L = L_iso(K_eff, mu_eff, "Kmu"); mat DM = M_M - M_A; //definition of the CTE tensor diff --git a/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp b/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp index b5274a769..d4172b5d1 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/SMA/unified_T.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -39,8 +40,9 @@ using namespace arma; namespace simcoon{ -void umat_sma_unified_T(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - +void umat_sma_unified_T(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -179,7 +181,7 @@ void umat_sma_unified_T(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat M_A = M_iso(K_A, mu_A, "Kmu"); mat M_M = M_iso(K_M, mu_M, "Kmu"); mat M = M_iso(K_eff, mu_eff, "Kmu"); - mat L = L_iso(K_eff, mu_eff, "Kmu"); + L = L_iso(K_eff, mu_eff, "Kmu"); mat DM = M_M - M_A; //definition of the CTE tensor diff --git a/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.cpp b/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.cpp index 488ced93c..36e7e1b71 100644 --- a/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Prony_Nfast.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -30,9 +31,10 @@ using namespace arma; namespace simcoon { -void umat_prony_Nfast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +void umat_prony_Nfast(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -62,6 +64,7 @@ void umat_prony_Nfast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, co //Define the viscoelastic stiffness mat L0 = L_iso(E0, nu0, "Enu"); + L = L0; mat M0 = M_iso(E0, nu0, "Enu"); ///@brief Temperature initialization double T_init = statev(0); diff --git a/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_Nfast.cpp b/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_Nfast.cpp index b26b0762b..acd389487 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_Nfast.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_Nfast.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -33,9 +34,10 @@ using namespace arma; namespace simcoon { -void umat_zener_Nfast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +void umat_zener_Nfast(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -65,6 +67,7 @@ void umat_zener_Nfast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, co //Define the viscoelastic stiffness mat L0 = L_iso(E0, nu0, "Enu"); + L = L0; ///@brief Temperature initialization double T_init = statev(0); diff --git a/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_fast.cpp b/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_fast.cpp index 462fd6d3b..6b6b23d48 100755 --- a/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_fast.cpp +++ b/src/Continuum_mechanics/Umat/Mechanical/Viscoelasticity/Zener_fast.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -30,9 +31,10 @@ using namespace arma; namespace simcoon { -void umat_zener_fast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) +void umat_zener_fast(const string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { - + + UNUSED(umat_name); UNUSED(nprops); UNUSED(nstatev); UNUSED(Time); @@ -72,6 +74,7 @@ void umat_zener_fast(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, con //Define the viscoelastic stiffness mat L0 = L_iso(E0, nu0, "Enu"); + L = L0; mat L1 = L_iso(E1, nu1, "Enu"); mat H1 = H_iso(etaB1, etaS1); //dimension of stiffness tensor diff --git a/src/Continuum_mechanics/Umat/umat_smart.cpp b/src/Continuum_mechanics/Umat/umat_smart.cpp index 1a5974c0b..32b9dc541 100644 --- a/src/Continuum_mechanics/Umat/umat_smart.cpp +++ b/src/Continuum_mechanics/Umat/umat_smart.cpp @@ -59,7 +59,6 @@ #include #include #include -#include #include #include #include @@ -263,35 +262,35 @@ void select_umat_M_finite(phase_characteristics &rve, const mat &DR,const double break; } case 2: { - umat_elasticity_iso(umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_elasticity_iso(rve.sptr_matprops->umat_name, umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 3: { - umat_elasticity_trans_iso(umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_elasticity_trans_iso(rve.sptr_matprops->umat_name, umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 4: { - umat_elasticity_ortho(umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_elasticity_ortho(rve.sptr_matprops->umat_name, umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 5: { - umat_hypoelasticity_ortho(umat_M->etot, umat_M->Detot, umat_M->F0, umat_M->F1, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_hypoelasticity_ortho(rve.sptr_matprops->umat_name, umat_M->etot, umat_M->Detot, umat_M->F0, umat_M->F1, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; - } + } case 6: { - umat_plasticity_iso_CCP(umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_plasticity_iso_CCP(rve.sptr_matprops->umat_name, umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 7: { - umat_plasticity_kin_iso_CCP(umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_plasticity_kin_iso_CCP(rve.sptr_matprops->umat_name, umat_M->etot, umat_M->Detot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 8: { - umat_saint_venant(umat_M->etot, umat_M->Detot, umat_M->F0, umat_M->F1, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_saint_venant(rve.sptr_matprops->umat_name, umat_M->etot, umat_M->Detot, umat_M->F0, umat_M->F1, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 9: { - umat_neo_hookean_incomp(umat_M->etot, umat_M->Detot, umat_M->F0, umat_M->F1, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + umat_neo_hookean_incomp(rve.sptr_matprops->umat_name, umat_M->etot, umat_M->Detot, umat_M->F0, umat_M->F1, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 10: case 11: case 12: case 13: case 14: case 15: { @@ -322,7 +321,7 @@ void select_umat_M(phase_characteristics &rve, const mat &DR,const double &Time, switch (list_umat[rve.sptr_matprops->umat_name]) { case 0: { - //umat_external(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + //umat_external(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); static dylib::library ext_lib("external/umat_plugin_ext", dylib::decorations::os_default()); @@ -337,7 +336,7 @@ void select_umat_M(phase_characteristics &rve, const mat &DR,const double &Time, ext_destroy ); - external_umat->umat_external_M(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + external_umat->umat_external_M(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } @@ -360,83 +359,80 @@ void select_umat_M(phase_characteristics &rve, const mat &DR,const double &Time, break; } case 2: { - umat_elasticity_iso(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_elasticity_iso(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 3: { - umat_elasticity_trans_iso(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_elasticity_trans_iso(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 4: { - umat_elasticity_ortho(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_elasticity_ortho(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 5: { - umat_plasticity_iso_CCP(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_plasticity_iso_CCP(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 6: { - umat_plasticity_kin_iso_CCP(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_plasticity_kin_iso_CCP(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 7: { - umat_plasticity_chaboche_CCP(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_plasticity_chaboche_CCP(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 8: { - umat_sma_unified_T(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + umat_sma_unified_T(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 9: { - umat_sma_aniso_T(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + umat_sma_aniso_T(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; - } + } case 10: { - umat_damage_LLD_0(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_damage_LLD_0(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 11: { - umat_zener_fast(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + umat_zener_fast(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 12: { - umat_zener_Nfast(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + umat_zener_Nfast(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 13: { - umat_prony_Nfast(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + umat_prony_Nfast(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 17: { - umat_plasticity_hill_isoh_CCP(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + umat_plasticity_hill_isoh_CCP(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 18: { - umat_hill_chaboche_CCP(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_hill_chaboche_CCP(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 19: { - umat_ani_chaboche_CCP(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_ani_chaboche_CCP(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 20: { - umat_dfa_chaboche_CCP(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_dfa_chaboche_CCP(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 21: { - umat_generic_chaboche_CCP(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, umat_M->sigma_in, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, solver_type, tnew_dt); + umat_generic_chaboche_CCP(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; - } + } case 22: { - umat_plasticity_hill_isoh_CCP_N(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + umat_plasticity_hill_isoh_CCP_N(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } - case 23: { - umat_sma_mono(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); - break; - } - case 24: { - umat_sma_mono_cubic(umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); + case 23: case 24: { + // SMAMO (isotropic) and SMAMC (cubic) both use unified umat_sma_mono + umat_sma_mono(rve.sptr_matprops->umat_name, umat_M->Etot, umat_M->DEtot, umat_M->sigma, umat_M->Lt, umat_M->L, DR, rve.sptr_matprops->nprops, rve.sptr_matprops->props, umat_M->nstatev, umat_M->statev, umat_M->T, umat_M->DT, Time, DTime, umat_M->Wm(0), umat_M->Wm(1), umat_M->Wm(2), umat_M->Wm(3), ndi, nshr, start, tnew_dt); break; } case 100: case 101: case 103: case 104: { diff --git a/testBin/Umats/UMEXT/external/umat_plugin_ext.cpp b/testBin/Umats/UMEXT/external/umat_plugin_ext.cpp index 84a00920a..573b2826c 100755 --- a/testBin/Umats/UMEXT/external/umat_plugin_ext.cpp +++ b/testBin/Umats/UMEXT/external/umat_plugin_ext.cpp @@ -26,8 +26,9 @@ class LIB_EXPORT umat_plugin_ext : public umat_plugin_ext_api { return "umext"; } - void umat_external_M(const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, vec &sigma_in, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, const int &solver_type, double &tnew_dt) { - + void umat_external_M(const std::string &umat_name, const vec &Etot, const vec &DEtot, vec &sigma, mat &Lt, mat &L, const mat &DR, const int &nprops, const vec &props, const int &nstatev, vec &statev, const double &T, const double &DT, const double &Time, const double &DTime, double &Wm, double &Wm_r, double &Wm_ir, double &Wm_d, const int &ndi, const int &nshr, const bool &start, double &tnew_dt) { + + UNUSED(umat_name); UNUSED(Etot); UNUSED(DR); UNUSED(nprops); @@ -37,9 +38,9 @@ class LIB_EXPORT umat_plugin_ext : public umat_plugin_ext_api { UNUSED(DTime); UNUSED(nshr); UNUSED(tnew_dt); - + double T_init = statev(0); - + //From the props to the material properties double E = props(0); double nu = props(1); @@ -47,38 +48,33 @@ class LIB_EXPORT umat_plugin_ext : public umat_plugin_ext_api { //Elastic stiffness tensor L = simcoon::L_iso(E, nu, "Enu"); - + ///@brief Initialization if(start) { T_init = T; sigma = zeros(6); - + Wm = 0.; Wm_r = 0.; Wm_ir = 0.; Wm_d = 0.; } - + vec sigma_start = sigma; - - //Compute the elastic strain and the related stress + + //Compute the elastic strain and the related stress vec Eel = Etot + DEtot - alpha*(T+DT-T_init); sigma = simcoon::el_pred(L, Eel, ndi); - - if((solver_type == 0)||(solver_type==2)) { - Lt = L; - } - else if(solver_type == 1) { - sigma_in = zeros(6); - } - + + Lt = L; + //Computation of the mechanical and thermal work quantities Wm += 0.5*sum((sigma_start+sigma)%DEtot); Wm_r += 0.5*sum((sigma_start+sigma)%DEtot); Wm_ir += 0.; Wm_d += 0.; - + statev(0) = T_init; }