diff --git a/.github/workflows/run-benchmark-examples.yml b/.github/workflows/run-benchmark-examples.yml deleted file mode 100644 index 45f9ec2..0000000 --- a/.github/workflows/run-benchmark-examples.yml +++ /dev/null @@ -1,65 +0,0 @@ -name: Examples-CI -on: - push: - - pull_request: - branches: [ main ] - - # Allows you to run this workflow manually from the Actions tab - workflow_dispatch: - - # Runs the workflow once per day at 3:15am - schedule: - - cron: '3 16 * * *' - -env: - CACHE_NUMBER: 1 # increase to reset cache manually - -jobs: - run-examples: - runs-on: ubuntu-latest - steps: - - name: checkout repo content - uses: actions/checkout@v2 - - - name: Setup Mambaforge - uses: conda-incubator/setup-miniconda@v3 - with: - miniforge-version: latest - activate-environment: model-validation - use-mamba: true - - - name: Set strict channel priority - run: conda config --set channel_priority strict - - - name: Update environment - run: mamba env update -n model-validation -f environment_benchmarks.yml - - - name: Run linear-elastic-plate-with-hole using fenics - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/examples/linear-elastic-plate-with-hole/fenics/ - python run_benchmark.py - - - name: Run linear-elastic-plate-with-hole using kratos - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/examples/linear-elastic-plate-with-hole/kratos/ - python run_benchmark.py - - - name: Archive results of the fenics run of linear-elastic-plate-with-hole - uses: actions/upload-artifact@v4 - with: - name: fenics_results_linear-elastic-plate-with-hole - path: | - examples/linear-elastic-plate-with-hole/fenics/results/ - - - name: Archive results of the kratos run of linear-elastic-plate-with-hole - uses: actions/upload-artifact@v4 - with: - name: kratos_results_linear-elastic-plate-with-hole - path: | - examples/linear-elastic-plate-with-hole/kratos/results/ - - - diff --git a/.github/workflows/run-benchmark.yml b/.github/workflows/run-benchmark.yml deleted file mode 100644 index 546edb1..0000000 --- a/.github/workflows/run-benchmark.yml +++ /dev/null @@ -1,160 +0,0 @@ -name: CI -on: - push: - - pull_request: - branches: [ main ] - - # Allows you to run this workflow manually from the Actions tab - workflow_dispatch: - - # Runs the workflow once per day at 3:15am - schedule: - - cron: '3 16 * * *' - -env: - CACHE_NUMBER: 1 # increase to reset cache manually - SNAKEMAKE_RESULT_FILE: metadata4ing_provenance - PROVENANACE_FILE_NAME: element_size_vs_max_mises_stress.pdf - -jobs: - run-simulation: - runs-on: ubuntu-latest - steps: - - name: checkout repo content - uses: actions/checkout@v2 - - - name: Setup Mambaforge - uses: conda-incubator/setup-miniconda@v3 - with: - miniforge-version: latest - activate-environment: model-validation - use-mamba: true - - - name: Set strict channel priority - run: conda config --set channel_priority strict - - - name: Update environment - run: mamba env update -n model-validation -f environment_benchmarks.yml - - - name: generate-config-files - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/benchmarks/linear-elastic-plate-with-hole/ - python generate_config.py - - - name: run_linear-elastic-plate-with-hole-benchmarks_snakemake - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/benchmarks/linear-elastic-plate-with-hole/ - snakemake --use-conda --force --cores 'all' - snakemake --use-conda --force --cores all \ - --reporter metadata4ing \ - --report-metadata4ing-paramscript ../common/parameter_extractor.py \ - --report-metadata4ing-config metadata4ing.config \ - --report-metadata4ing-filename $SNAKEMAKE_RESULT_FILE - - - name: run_linear-elastic-plate-with-hole-benchmarks_nextflow - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/benchmarks/linear-elastic-plate-with-hole/ - nextflow run main.nf -params-file workflow_config.json -c ../common/nextflow.config -plugins nf-prov@1.4.0 - - - name: Archive Linear Elastic plate with a hole benchmark data for snakemake - uses: actions/upload-artifact@v4 - with: - name: snakemake_results_linear-elastic-plate-with-hole - path: | - benchmarks/linear-elastic-plate-with-hole/${{ env.SNAKEMAKE_RESULT_FILE }}.zip - - - name: Archive Linear Elastic plate with a hole benchmark data for nextflow - uses: actions/upload-artifact@v4 - with: - name: nextflow_results_linear-elastic-plate-with-hole - path: | - benchmarks/linear-elastic-plate-with-hole/nextflow_results/ - - process-artifacts: - runs-on: ubuntu-latest - needs: run-simulation - steps: - - name: Checkout repo content - uses: actions/checkout@v2 - - - name: Download artifact - uses: actions/download-artifact@v4 - with: - name: snakemake_results_linear-elastic-plate-with-hole - path: ./artifact_files - - - name: Unzip Snakemake Result File - run: | - mkdir -p ./$SNAKEMAKE_RESULT_FILE - unzip -o ./artifact_files/$SNAKEMAKE_RESULT_FILE.zip -d ./$SNAKEMAKE_RESULT_FILE - - - name: Setup Mambaforge with postprocessing env - uses: conda-incubator/setup-miniconda@v3 - with: - miniforge-version: latest - activate-environment: postprocessing - use-mamba: true - environment-file: benchmarks/linear-elastic-plate-with-hole/environment_postprocessing.yml - - - name: Validate Snakemake Result File - shell: bash -l {0} - run: | - python benchmarks/common/validate_provenance.py \ - --provenance_folderpath "./$SNAKEMAKE_RESULT_FILE" - - - name: Run plotting script - shell: bash -l {0} - run: | - python benchmarks/linear-elastic-plate-with-hole/plot_metrics.py \ - --provenance_folderpath "./$SNAKEMAKE_RESULT_FILE" \ - --output_file $PROVENANACE_FILE_NAME - - - name: Upload snakemake results file as artifact - uses: actions/upload-artifact@v4 - with: - name: element-size-vs-stress-plot - path: ${{ env.PROVENANACE_FILE_NAME }} - - - name: Re-zip snakemake result folder - run: | - cd "./${SNAKEMAKE_RESULT_FILE}" - zip -r "../${SNAKEMAKE_RESULT_FILE}.zip" . - - - name: Upload RoCrate Zip file onto RoHub - id: rohub_upload - shell: bash -l {0} - continue-on-error: true - run: | - LOG_FILE="${GITHUB_WORKSPACE}/rohub_upload_error.log" - - # Run Python and capture stdout+stderr - OUTPUT="$( - python benchmarks/common/upload_provenance.py \ - --provenance_folderpath "./${SNAKEMAKE_RESULT_FILE}.zip" \ - --benchmark_name "linear-elastic-plate-with-hole" \ - --username "${{ secrets.ROHUB_USERNAME }}" \ - --password "${{ secrets.ROHUB_PASSWORD }}" \ - 2>&1 - )" - PYTHON_EXIT_CODE=$? - - # Standard convention: 0 = success, non-zero = failure - if [ "$PYTHON_EXIT_CODE" -eq 0 ]; then - echo "RoHub upload succeeded" - rm -f "$LOG_FILE" || true - else - echo "=== RoHub upload failed — writing log to $LOG_FILE ===" - printf '%s\n' "$OUTPUT" | tee "$LOG_FILE" - fi - - - name: Upload RoHub error log - if: always() - uses: actions/upload-artifact@v4 - with: - name: rohub-upload-error-log - path: ${{ github.workspace }}/rohub_upload_error.log - if-no-files-found: ignore diff --git a/.github/workflows/run-dumux-benchmark.yml b/.github/workflows/run-dumux-benchmark.yml deleted file mode 100644 index 04dddf5..0000000 --- a/.github/workflows/run-dumux-benchmark.yml +++ /dev/null @@ -1,51 +0,0 @@ -name: DuMux-CI -on: - push: - - pull_request: - branches: [ main ] - - workflow_dispatch: - -env: - CACHE_NUMBER: 1 # increase to reset cache manually - -jobs: - tests: - runs-on: ubuntu-latest - - steps: - - name: checkout repo content - uses: actions/checkout@v2 - - - name: Setup Mambaforge - uses: conda-incubator/setup-miniconda@v3 - with: - miniforge-version: latest - activate-environment: model-validation - use-mamba: true - - - name: Set strict channel priority - run: conda config --set channel_priority strict - - - name: Setup Apptainer - uses: eWaterCycle/setup-apptainer@v2 - with: - apptainer-version: 1.4.5 - - - name: Update environment - run: mamba env update -n model-validation -f environment_benchmarks.yml - - - name: Simulate rotating cylinders using dumux & validate convergence - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/benchmarks/rotating-cylinders/dumux - python3 run_benchmark.py - python3 convergence_test.py - - - name: Upload results - uses: actions/upload-artifact@v4 - with: - name: convergence-results - path: benchmarks/rotating-cylinders/dumux/rotating-cylinders/results - retention-days: 5 \ No newline at end of file diff --git a/.github/workflows/run-julia-benchmark.yml b/.github/workflows/run-edelweiss-benchmark.yml similarity index 65% rename from .github/workflows/run-julia-benchmark.yml rename to .github/workflows/run-edelweiss-benchmark.yml index 748a39b..b5d2833 100644 --- a/.github/workflows/run-julia-benchmark.yml +++ b/.github/workflows/run-edelweiss-benchmark.yml @@ -31,24 +31,27 @@ jobs: miniforge-version: latest activate-environment: model-validation use-mamba: true + - name: Setup Apptainer uses: eWaterCycle/setup-apptainer@v2 with: apptainer-version: 1.4.5 + - name: Set strict channel priority run: conda config --set channel_priority strict - name: Update environment run: mamba env update -n model-validation -f environment_benchmarks.yml - - name: generate-config-files + - name: Run linear-elastic-plate-with-hole using edelweiss shell: bash -l {0} run: | - cd $GITHUB_WORKSPACE/benchmarks/linear-elastic-plate-with-hole/ - python generate_config.py + cd $GITHUB_WORKSPACE/examples/linear-elastic-plate-with-hole/edelweiss/ + snakemake --use-conda --force --cores 'all' --use-apptainer - - name: run_linear-elastic-plate-with-hole-benchmarks_snakemake - shell: bash -l {0} - run: | - cd $GITHUB_WORKSPACE/benchmarks/linear-elastic-plate-with-hole/ - snakemake --use-conda --use-apptainer --cores=all --config tools=[extendablefem] configurations=["1","05","025","0125","00625","003125"] \ No newline at end of file + - name: Archive results of the edelweiss run of linear-elastic-plate-with-hole + uses: actions/upload-artifact@v4 + with: + name: edelweiss_results_linear-elastic-plate-with-hole + path: | + examples/linear-elastic-plate-with-hole/edelweiss/solution_metrics.json \ No newline at end of file diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/README.md b/examples/linear-elastic-plate-with-hole/edelweiss/README.md new file mode 100755 index 0000000..c0d2445 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/README.md @@ -0,0 +1,8 @@ +# EdelweissFE_Snakemake-Workflow + +## Getting started + +``` +mamba install -c bioconda snakemake +snakemake --use-conda --cores all +``` diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/Snakefile b/examples/linear-elastic-plate-with-hole/edelweiss/Snakefile new file mode 100755 index 0000000..e010073 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/Snakefile @@ -0,0 +1,74 @@ +import json + +rule all: + input: + "esExport.case", + "solution_metrics.json" + +rule get_EdelweissFE: + params: + target = "https://github.com/Edelweiss-Numerics/EdelweissFE-Apptainer/releases/download/v26.03/edelweissfe-debian-13.sif", + output: + container = "edelweiss.sif", + shell: + """ + curl -L -o {output.container} {params.target} + """ + +rule create_mesh: + input: + parameters = "parameters.json", + output: + mesh_temp = "mesh.msh", + mesh = "mesh.inp", + conda: "environment_mesh.yml" + shell: + """ + python3 "create_mesh.py" --input_parameter_file {input.parameters} --output_mesh_file {output.mesh_temp} + # create {output.mesh} from {output.mesh_temp} + python3 convert_mesh.py + # meshio c {output.mesh_temp} {output.mesh} + # remove heading + sed -i '1,3d' {output.mesh} + # change element type to CPE4 + sed -i 's/CAX4P/CPE4/' {output.mesh} + """ + +rule prepare_simulation: + input: + "parameters.json", + "analytical_solution.py", + "input_template.inp" + output: + "input.inp", + "boundaryFields.inp" + conda: "environment_preProcessing.yml" + shell: + """ + python3 create_BCsFromAnalyticalSolution.py + python3 create_inputFromInputTemplate.py + """ + +rule run_edelweiss: + input: + "edelweiss.sif", + "input.inp", + "mesh.inp", + "boundaryFields.inp" + output: + "esExport.case" + shell: + """ + apptainer exec --bind $(pwd) edelweiss.sif edelweissfe input.inp + """ + +rule get_results: + input: + ensight = "esExport.case" + output: + metrics = "solution_metrics.json", + conda: "environment_postProcessing.yml" + shell: + """ + python3 create_metrics.py {input.ensight} {output.metrics} + """ diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/analytical_solution.py b/examples/linear-elastic-plate-with-hole/edelweiss/analytical_solution.py new file mode 100755 index 0000000..6ff7fe4 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/analytical_solution.py @@ -0,0 +1,144 @@ +import numpy as np +import sympy as sp + +class AnalyticalSolution: + def __init__(self, E: float, nu: float, radius: float, L:float, load:float) -> None: + self.radius = radius + self.L = L + self.load = load + self.E = E + self.nu = nu + + # Use symbolic expressions from displacement_symbolic + # to create fast lambdified functions for numerical evaluation + ux_sym, uy_sym = self.displacement_symbolic() + x, y = sp.symbols('x y') + E, nu, a, T = sp.symbols('E nu a T') + subs = { + E: self.E, + nu: self.nu, + a: self.radius, + T: self.load + } + ux_expr_num = ux_sym.subs(subs) + uy_expr_num = uy_sym.subs(subs) + self.ux_func = sp.lambdify((x, y), ux_expr_num, modules='numpy') + self.uy_func = sp.lambdify((x, y), uy_expr_num, modules='numpy') + + # Symbolic stress expressions and lambdify + sxx_sym, sxy_sym, syy_sym = self.stress_symbolic() + x, y = sp.symbols('x y') + E, nu, a, T = sp.symbols('E nu a T') + subs = { + E: self.E, + nu: self.nu, + a: self.radius, + T: self.load + } + sxx_expr_num = sxx_sym.subs(subs) + sxy_expr_num = sxy_sym.subs(subs) + syy_expr_num = syy_sym.subs(subs) + self.sxx_func = sp.lambdify((x, y), sxx_expr_num, modules='numpy') + self.sxy_func = sp.lambdify((x, y), sxy_expr_num, modules='numpy') + self.syy_func = sp.lambdify((x, y), syy_expr_num, modules='numpy') + + def displacement_symbolic(self): + """ + Returns symbolic expressions for ux and uy in terms of x and y. + """ + # Define symbols + x, y = sp.symbols('x y') + E, nu, a, T = sp.symbols('E nu a T') + r = sp.sqrt(x**2 + y**2) + theta = sp.atan2(y, x) + Ta_8mu = T * a / (4 * E / (1.0 + nu)) + k = (3.0 - nu) / (1.0 + nu) + + ct = sp.cos(theta) + c3t = sp.cos(3 * theta) + st = sp.sin(theta) + s3t = sp.sin(3 * theta) + + fac = 2 * (a / r)**3 + + ux = Ta_8mu * ( + r / a * (k + 1.0) * ct + 2.0 * a / r * ((1.0 + k) * ct + c3t) - fac * c3t + ) + + uy = Ta_8mu * ( + (r / a) * (k - 3.0) * st + 2.0 * a / r * ((1.0 - k) * st + s3t) - fac * s3t + ) + + return ux, uy + + def stress_symbolic(self): + """ + Returns symbolic expressions for sxx, sxy, syx, syy in terms of x and y. + """ + x, y = sp.symbols('x y') + E, nu, a, T = sp.symbols('E nu a T') + r = sp.sqrt(x**2 + y**2) + theta = sp.atan2(y, x) + cos2t = sp.cos(2 * theta) + cos4t = sp.cos(4 * theta) + sin2t = sp.sin(2 * theta) + sin4t = sp.sin(4 * theta) + + fac1 = (a * a) / (r * r) + fac2 = 1.5 * fac1 * fac1 + + sxx = T - T * fac1 * (1.5 * cos2t + cos4t) + T * fac2 * cos4t + syy = -T * fac1 * (0.5 * cos2t - cos4t) - T * fac2 * cos4t + sxy = -T * fac1 * (0.5 * sin2t + sin4t) + T * fac2 * sin4t + + return sxx, sxy, syy + + def displacement(self, x: np.ndarray) -> np.ndarray: + """ + Evaluates the symbolic displacement expressions + Accepts x of shape (2, N) or (3,N) (third dimension is omitted, but some tools always compute in 3D) + and returns two 1D arrays. + The lambda functions are computed from the symbolic representation in the constructor + """ + arr = np.asarray(x) + if arr.ndim != 2 or arr.shape[0] not in (2, 3): + raise ValueError(f"Input x must have shape (2, N) or (3, N), got {arr.shape}") + # If 3D, ignore the third row + ux_vals = self.ux_func(arr[0], arr[1]) + uy_vals = self.uy_func(arr[0], arr[1]) + result = ux_vals, uy_vals + + return result + + def displacement_symbolic_str(self, X_str: str, Y_str: str): + """ + Returns string representations of the symbolic displacement functions + with variable names x and y being replaced by X_str and Y_str. + """ + # Get symbolic expressions + ux, uy = self.displacement_symbolic() + # Define symbols + x, y = sp.symbols('x y') + E, nu, a, T = sp.symbols('E nu a T') + subs_vars = {x: sp.Symbol(X_str), y: sp.Symbol(Y_str)} + subs_params = {E: self.E, nu: self.nu, a: self.radius, T: self.load} + ux_sub = ux.subs(subs_vars).subs(subs_params) + uy_sub = uy.subs(subs_vars).subs(subs_params) + # Convert to string using sympy.sstr for single-line output with ** + ux_str = sp.sstr(ux_sub) + uy_str = sp.sstr(uy_sub) + return ux_str, uy_str + + def stress(self, x: np.ndarray) -> np.ndarray: + """ + Evaluates the symbolic stress expressions. + Accepts coordinates x of shape (2, N) or (3, N) (third dimension is omitted). + Returns three 1D arrays: sxx, sxy, syy. + """ + arr = np.asarray(x) + if arr.ndim != 2 or arr.shape[0] not in (2, 3): + raise ValueError(f"Input x must have shape (2, N) or (3, N), got {arr.shape}") + sxx = self.sxx_func(arr[0], arr[1]) + sxy = self.sxy_func(arr[0], arr[1]) + syy = self.syy_func(arr[0], arr[1]) + return sxx, sxy, syy diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/convert_mesh.py b/examples/linear-elastic-plate-with-hole/edelweiss/convert_mesh.py new file mode 100755 index 0000000..2e74308 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/convert_mesh.py @@ -0,0 +1,62 @@ +import meshio +import numpy as np + +mesh = meshio.read("mesh.msh") + +# Keep only 2D cells +cells = [] +cell_data = { + "abaqus:element_type": [] +} + +for cell_block in mesh.cells: + if cell_block.type == "triangle": + cells.append(cell_block) + cell_data["abaqus:element_type"].append( + ["CPS3"] * len(cell_block.data) + ) + + elif cell_block.type == "quad": + cells.append(cell_block) + cell_data["abaqus:element_type"].append( + ["CPS4"] * len(cell_block.data) + ) + +# Use first and second coordinates only +pts2d = mesh.points[:, :2] +c1 = pts2d[:, 0] +c2 = pts2d[:, 1] + +# Bounds +c1min, c1max = float(c1.min()), float(c1.max()) +c2min, c2max = float(c2.min()), float(c2.max()) + +# Tolerances (relative to span, with a tiny absolute floor) +c1tol = max(1e-12, 1e-8 * max(1.0, c1max - c1min)) +c2tol = max(1e-12, 1e-8 * max(1.0, c2max - c2min)) + +# Indices for node sets +i_c1min = np.where(np.isclose(c1, c1min, atol=c1tol))[0] +i_c1max = np.where(np.isclose(c1, c1max, atol=c1tol))[0] +i_c2min = np.where(np.isclose(c2, c2min, atol=c2tol))[0] +i_c2max = np.where(np.isclose(c2, c2max, atol=c2tol))[0] + +point_sets = { + "X_MIN": i_c1min.astype(np.int32), + "X_MAX": i_c1max.astype(np.int32), + "Y_MIN": i_c2min.astype(np.int32), + "Y_MAX": i_c2max.astype(np.int32), +} + +mesh2d = meshio.Mesh( + points=pts2d, + cells=cells, + cell_data=cell_data, + point_sets=point_sets, +) + + + +meshio.write("mesh.inp", mesh2d) + + diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/create_BCsFromAnalyticalSolution.py b/examples/linear-elastic-plate-with-hole/edelweiss/create_BCsFromAnalyticalSolution.py new file mode 100755 index 0000000..0174ba2 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/create_BCsFromAnalyticalSolution.py @@ -0,0 +1,41 @@ +import json + +from analytical_solution import AnalyticalSolution + +with open("parameters.json", "r", encoding="utf-8") as file: + parameters = json.load(file) + +assert parameters["young-modulus"]["unit"] == "Pa" +assert parameters["poisson-ratio"]["unit"] == "" +assert parameters["radius"]["unit"] == "m" +assert parameters["length"]["unit"] == "m" +assert parameters["load"]["unit"] == "MPa" + +E = parameters["young-modulus"]["value"] +nu = parameters["poisson-ratio"]["value"] +R = parameters["radius"]["value"] +L = parameters["length"]["value"] +load = parameters["load"]["value"] * 1e6 # convert to Pa + +solutionSymbolic = AnalyticalSolution(E, nu, R, L, load) + +def replaceMathFunctions(expression): + expression = expression.replace("sqrt", "np.sqrt") + expression = expression.replace("sin", "np.sin") + expression = expression.replace("cos", "np.cos") + expression = expression.replace("atan2", "np.atan2") + return expression + +with open("boundaryFields.inp", "w+") as f: + f.write("*analyticalField, name=ux, type=scalarExpression" + "\n") + expression = solutionSymbolic.displacement_symbolic_str("x", "y")[0] + expression = replaceMathFunctions(expression) + f.write( + '"f(x,y,z)"="' + expression + '"\n' + ) + f.write("*analyticalField, name=uy, type=scalarExpression" + "\n") + expression = solutionSymbolic.displacement_symbolic_str("x", "y")[1] + expression = replaceMathFunctions(expression) + f.write( + '"f(x,y,z)"="' + expression + '"\n' + ) diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/create_inputFromInputTemplate.py b/examples/linear-elastic-plate-with-hole/edelweiss/create_inputFromInputTemplate.py new file mode 100755 index 0000000..60cd2a9 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/create_inputFromInputTemplate.py @@ -0,0 +1,19 @@ +import json + +with open("parameters.json", "r", encoding="utf-8") as file: + parameters = json.load(file) + +with open("./input_template.inp", "r") as infile: + lines = infile.readlines() + + +replaceDict = { + "_E_": parameters["young-modulus"]["value"], + "_NU_": parameters["poisson-ratio"]["value"], +} + +with open("input.inp", "w+") as outfile: + for line in lines: + for old, new in replaceDict.items(): + line = line.replace(old,str(new)) + outfile.write(line) diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/create_mesh.py b/examples/linear-elastic-plate-with-hole/edelweiss/create_mesh.py new file mode 100755 index 0000000..5a93d81 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/create_mesh.py @@ -0,0 +1,106 @@ +import json +from argparse import ArgumentParser +import os + +import gmsh +from pint import UnitRegistry + +ureg = UnitRegistry() + + +def create_mesh(parameter_file, mesh_file): + # Load parameters + with open(parameter_file) as f: + parameters = json.load(f) + print(parameters) + + # Read configuration from parameters instead of the filename + configuration = parameters["configuration"] + + length = ( + ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"]) + .to_base_units() + .magnitude + ) + radius = ( + ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"]) + .to_base_units() + .magnitude + ) + # create mesh with gmsh python api + r""" + 4---------3 + | | + 5_ | + \ | + 1______2 + + """ + + gmsh.initialize() + gmsh.model.add(configuration) + + element_size = ( + ureg.Quantity( + parameters["element-size"]["value"], parameters["element-size"]["unit"] + ) + .to_base_units() + .magnitude + ) + + gmsh.option.setNumber("Mesh.CharacteristicLengthMin", element_size) + gmsh.option.setNumber("Mesh.CharacteristicLengthMax", element_size) + gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 1.0) + gmsh.option.setNumber("Mesh.ElementOrder", parameters["element-order"]) + + z = 0.0 + lc = 1.0 + + x0 = 0.0 + x1 = x0 + radius + x2 = x0 + length + y0 = 0.0 + y1 = y0 + radius + y2 = y0 + length + + center = gmsh.model.geo.addPoint(x0, y0, z, lc) + p1 = gmsh.model.geo.addPoint(x1, y0, z, lc) + p2 = gmsh.model.geo.addPoint(x2, y0, z, lc) + p3 = gmsh.model.geo.addPoint(x2, y2, z, lc) + p4 = gmsh.model.geo.addPoint(x0, y2, z, lc) + p5 = gmsh.model.geo.addPoint(x0, y1, z, lc) + + l1 = gmsh.model.geo.addLine(p1, p2) + l2 = gmsh.model.geo.addLine(p2, p3) + l3 = gmsh.model.geo.addLine(p3, p4) + l4 = gmsh.model.geo.addLine(p4, p5) + l5 = gmsh.model.geo.addCircleArc(p5, center, p1) + + curve = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4, l5]) + plane = gmsh.model.geo.addPlaneSurface([curve]) + gmsh.model.geo.synchronize() + gmsh.model.geo.removeAllDuplicates() + gmsh.model.addPhysicalGroup(2, [plane], 1, name="surface") + gmsh.model.addPhysicalGroup(1, [l4], 1, name="boundary_left") + gmsh.model.addPhysicalGroup(1, [l1], 2, name="boundary_bottom") + gmsh.model.addPhysicalGroup(1, [l2], 3, name="boundary_right") + gmsh.model.addPhysicalGroup(1, [l3], 4, name="boundary_top") + + # Recombine triangles into quads + gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) + gmsh.option.setNumber("Mesh.RecombineAll", 1) + gmsh.option.setNumber("Mesh.Algorithm", 8) + + gmsh.model.mesh.generate(2) + + + gmsh.write(mesh_file) + gmsh.finalize() + + +if __name__ == "__main__": + PARSER = ArgumentParser(description="Create input files and mesh for FEniCS simulation") + PARSER.add_argument("--input_parameter_file", required=True, help="JSON file containing simulation parameters") + PARSER.add_argument("--output_mesh_file", required=True, help="Output path for the generated mesh (.msh)") + ARGS = vars(PARSER.parse_args()) + create_mesh(ARGS["input_parameter_file"], ARGS["output_mesh_file"]) diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/create_metrics.py b/examples/linear-elastic-plate-with-hole/edelweiss/create_metrics.py new file mode 100755 index 0000000..c5cca63 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/create_metrics.py @@ -0,0 +1,47 @@ +import json +import sys + +import numpy as np +import pyvista as pv + +with open("parameters.json", "r", encoding="utf-8") as file: + parameters = json.load(file) + +caseFile = sys.argv[1] +outFile = sys.argv[2] + +reader = pv.EnSightReader(caseFile) +reader.set_active_time_set(1) +reader.set_active_time_point(len(reader.time_values) - 1) + +mesh = reader.read()["all"] + +voigtStress = mesh["stress"] +stressMises = np.sqrt( + voigtStress[:, 0] ** 2 + + voigtStress[:, 1] ** 2 + + voigtStress[:, 2] ** 2 + - voigtStress[:, 1] * voigtStress[:, 2] + - voigtStress[:, 0] * voigtStress[:, 2] + - voigtStress[:, 0] * voigtStress[:, 1] + + 3 * (voigtStress[:, 3] ** 2 + voigtStress[:, 4] ** 2 + voigtStress[:, 5] ** 2) +) + +x = parameters["displacement-evaluation-point"]["x"]["value"] +y = parameters["displacement-evaluation-point"]["y"]["value"] + +point = pv.PointSet([x, y, 0]) + +sampled = point.sample(mesh) + +metrics = { + # "max_von_mises_stress_nodes": "", + "max_von_mises_stress_gauss_points": float(stressMises.max()), + f"displacement_x_at_evaluation_point (x=f'{x:.2f}', y=f'{y:.2f}')": float( + sampled["displacement"][0][0] + ), +} + +with open(outFile, "w+") as f: + f.write(json.dumps(metrics, indent=4)) + # f.write() diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/environment_mesh.yml b/examples/linear-elastic-plate-with-hole/edelweiss/environment_mesh.yml new file mode 100755 index 0000000..4ebb03c --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/environment_mesh.yml @@ -0,0 +1,13 @@ +name: mesh-generation +# Environment file for create_mesh.py script. Called by the main workflow. + +channels: + - conda-forge + +channel_priority: strict + +dependencies: + - python=3.12 + - pint + - python-gmsh + - meshio diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/environment_postProcessing.yml b/examples/linear-elastic-plate-with-hole/edelweiss/environment_postProcessing.yml new file mode 100755 index 0000000..56a3dde --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/environment_postProcessing.yml @@ -0,0 +1,11 @@ +name: mesh-generation +# Environment file for create_mesh.py script. Called by the main workflow. + +channels: + - conda-forge + +channel_priority: strict + +dependencies: + - python=3.12 + - pyvista diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/environment_preProcessing.yml b/examples/linear-elastic-plate-with-hole/edelweiss/environment_preProcessing.yml new file mode 100755 index 0000000..6bc76b7 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/environment_preProcessing.yml @@ -0,0 +1,12 @@ +name: mesh-generation +# Environment file for create_mesh.py script. Called by the main workflow. + +channels: + - conda-forge + +channel_priority: strict + +dependencies: + - python=3.12 + - numpy + - sympy diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/input_template.inp b/examples/linear-elastic-plate-with-hole/edelweiss/input_template.inp new file mode 100755 index 0000000..42c9237 --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/input_template.inp @@ -0,0 +1,38 @@ +*job, name=test, domain=2d + +*solver, solver=NISTParallel, name=theSolver + +*material, name=LINEARELASTIC, id=mat +_E_, _NU_ + +*include, input="mesh.inp" + +*section, name=section1, material=mat, type=plane, thickness=1. +all + +*fieldOutput +create=perNode, name=displacement, elSet=all, field=displacement, result=U +create=perElement, name=stress, elSet=all, result=stress, quadraturePoint=0:4, f(x)='np.max(x,axis=1)' + +** create=perNode, name=RF, nSet=x_min, field=displacement, result=P, f(x)='sum(x[:,1])', saveHistory=True, export=RF + +***output, type=meshplot, name=RFPlot +**figure=1, axSpec=111, c=k, ls=solid, create=xyData, x=time, y=RF + +*output, type=ensight, name=esExport +create=perNode, fieldOutput=displacement +create=perElement, fieldOutput=stress +configuration, overwrite=yes + +*include, input="boundaryFields.inp" + +*step, maxInc=1e-0, minInc=1e-3, maxNumInc=1000, maxIter=25, stepLength=1, solver=theSolver +dirichlet, name=xmin, nSet=X_MIN, field=displacement, 1=0. +dirichlet, name=ymin, nSet=Y_MIN, field=displacement, 2=0. + +*step, maxInc=1e-0, minInc=1e-3, maxNumInc=1000, maxIter=25, stepLength=1, solver=theSolver +dirichlet, name=xmax_1, nSet=X_MAX, field=displacement, 1=1, analyticalField=ux +dirichlet, name=xmax_2, nSet=X_MAX, field=displacement, 2=1, analyticalField=uy +dirichlet, name=ymax_1, nSet=Y_MAX, field=displacement, 1=1, analyticalField=ux +dirichlet, name=ymax_2, nSet=Y_MAX, field=displacement, 2=1, analyticalField=uy + diff --git a/examples/linear-elastic-plate-with-hole/edelweiss/parameters.json b/examples/linear-elastic-plate-with-hole/edelweiss/parameters.json new file mode 100755 index 0000000..8b2ea2c --- /dev/null +++ b/examples/linear-elastic-plate-with-hole/edelweiss/parameters.json @@ -0,0 +1,41 @@ +{ + "configuration": "1", + "radius": { + "value": 0.33, + "unit": "m" + }, + "length": { + "value": 1.0, + "unit": "m" + }, + "displacement-evaluation-point": { + "x": { + "value": 0.33, + "unit": "m" + }, + "y": { + "value": 0.0, + "unit": "m" + } + }, + "load": { + "value": 100.0, + "unit": "MPa" + }, + "element-size": { + "value": 0.1, + "unit": "m" + }, + "element-order": 1, + "element-degree": 1, + "quadrature-rule": "gauss", + "quadrature-degree": 1, + "young-modulus": { + "value": 210e9, + "unit": "Pa" + }, + "poisson-ratio": { + "value": 0.3, + "unit": "" + } +} \ No newline at end of file