diff --git a/.github/workflows/arch.yml b/.github/workflows/arch.yml new file mode 100644 index 0000000..db412b8 --- /dev/null +++ b/.github/workflows/arch.yml @@ -0,0 +1,19 @@ +name: Arch +on: [pull_request] +jobs: + test: + runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + include: + - alternative: "" + label: Default + - alternative: "1" + label: Alternative + name: test (${{ matrix.label }}) + steps: + - uses: actions/checkout@v5 + - name: Build docker image + run: | + bash zscripts/docker-build-arch.bash ${{ matrix.alternative }} diff --git a/.github/workflows/test_and_coverage.yml b/.github/workflows/ubuntu.yml similarity index 72% rename from .github/workflows/test_and_coverage.yml rename to .github/workflows/ubuntu.yml index a4c1938..c5cb0b7 100644 --- a/.github/workflows/test_and_coverage.yml +++ b/.github/workflows/ubuntu.yml @@ -1,13 +1,13 @@ -name: Test & Coverage +name: Ubuntu on: [pull_request] jobs: - test_and_coverage: + test: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - name: Install Libraries run: | - sudo apt-get update -y && sudo apt-get install -y liblapacke-dev libmumps-seq-dev libopenblas-dev libsuitesparse-dev python3-matplotlib + sudo apt-get update -y && sudo apt-get install -y liblapacke-dev libopenblas-dev libsuitesparse-dev - name: Run tests run: | cargo test -- --nocapture @@ -16,7 +16,7 @@ jobs: curl -LsSf https://github.com/taiki-e/cargo-llvm-cov/releases/latest/download/cargo-llvm-cov-x86_64-unknown-linux-gnu.tar.gz | tar xzf - -C ~/.cargo/bin - name: Generate code coverage run: | - cargo llvm-cov --all-features --workspace --ignore-filename-regex build.rs --lcov --output-path lcov.info + cargo llvm-cov --workspace --ignore-filename-regex build.rs --lcov --output-path lcov.info - name: Upload to codecov.io uses: codecov/codecov-action@v3 with: diff --git a/.vscode/settings.json b/.vscode/settings.json index ff4734b..fc75949 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,19 +1,51 @@ { "cSpell.words": [ + "archlinux", "argsort", "Bbแตแตข", + "blas", "connectivities", "conv", + "cpmech", "cuthill", "deque", + "Dorival", "drawmsh", + "Durand", "extrap", + "Farias", + "gemlab", + "integ", + "ldconfig", + "liblapacke", + "libopenblas", + "libsuitesparse", + "makepkg", + "matplotlib", "nfacet", "ngauss", + "ngausss", "nmarked", + "nnode", + "nocapture", + "noconfirm", + "NOPASSWD", + "oneapi", + "openblas", + "pacman", + "Pedroso", "PINV", + "plotpy", + "serde", + "setvars", "Stockmeyer", + "structopt", + "suitesparse", + "tetgen", + "tritet", "Warshall", - "zrange" + "zdocker", + "zrange", + "zscripts" ] } \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml index 43f86f0..1787b6e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,16 +11,25 @@ readme = "README.md" categories = ["mathematics", "science"] keywords = ["geometry", "mesh", "fem"] +[features] +intel_mkl = [ + "russell_lab/intel_mkl", + "russell_sparse/intel_mkl", + "russell_stat/intel_mkl", + "russell_tensor/intel_mkl", +] +local_sparse = ["russell_sparse/local_sparse"] + [dependencies] -plotpy = "1.14" -russell_lab = "1.10" -russell_sparse = "1.10" -russell_stat = "1.10" -russell_tensor = "1.10" -tritet = "2.0" +plotpy = "1.23" +russell_lab = "2.1" +russell_sparse = "2.1" +russell_stat = "2.1" +russell_tensor = "2.1" +tritet = { version = "3.2", features = ["with_tetgen"] } serde = { version = "1.0", features = ["derive"] } serde_json = "1.0" structopt = "0.3" [dev-dependencies] -criterion = "0.7" +criterion = "0.8" diff --git a/README.md b/README.md index 544510a..10c2329 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,8 @@ # Geometry, meshes, and numerical integration for finite element analyses [![codecov](https://codecov.io/gh/cpmech/gemlab/graph/badge.svg?token=OjQKQ0PrNF)](https://codecov.io/gh/cpmech/gemlab) -[![Test & Coverage](https://github.com/cpmech/gemlab/actions/workflows/test_and_coverage.yml/badge.svg)](https://github.com/cpmech/gemlab/actions/workflows/test_and_coverage.yml) +[![Arch](https://github.com/cpmech/gemlab/actions/workflows/arch.yml/badge.svg)](https://github.com/cpmech/gemlab/actions/workflows/arch.yml) +[![Ubuntu](https://github.com/cpmech/gemlab/actions/workflows/ubuntu.yml/badge.svg)](https://github.com/cpmech/gemlab/actions/workflows/ubuntu.yml) [![documentation](https://img.shields.io/badge/gemlab-documentation-blue)](https://docs.rs/gemlab) @@ -10,9 +11,12 @@ - [Introduction](#introduction) - [Documentation](#documentation) - [Installation](#installation) - - [TL;DR (Debian/Ubuntu/Linux)](#tldr-debianubuntulinux) - - [Details](#details) - - [Setting Cargo.toml](#setting-cargotoml) + - [Arch Linux](#arch-linux) + - [1. Default](#1-default) + - [2. Alternative](#2-alternative) + - [Debian/Ubuntu Linux](#debianubuntu-linux) + - [1. Default](#1-default-1) + - [2. Alternative](#2-alternative-1) - [Examples](#examples) - [MSH file format](#msh-file-format) - [Numerical integration](#numerical-integration) @@ -38,44 +42,44 @@ This crate contains structures and functions for geometry computations, generate ## Installation -At this moment, Gemlab works on **Linux** (Debian/Ubuntu; and maybe Arch). +This crate depends on [Russell (Rust Scientific Library)](https://github.com/cpmech/russell) and, therefore, has the same dependencies as `russell`. Among the combinations described in `russell` website, the following are recommended here: -### TL;DR (Debian/Ubuntu/Linux) +1. **Default** - Use [OpenBLAS](https://github.com/OpenMathLib/OpenBLAS) and [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse) +2. **Alternative** - Use [Intel MKL](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) and locally compiled [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse) and the [MUMPS solver](https://mumps-solver.org) -First: +Therefore, the following (re-exported) **features** are available here: -```bash -sudo apt-get install -y --no-install-recommends \ - g++ \ - gdb \ - gfortran \ - liblapacke-dev \ - libmumps-seq-dev \ - libopenblas-dev \ - libsuitesparse-dev -``` +* `intel_mkl` - Tells `russell` to use Intel MKL +* `local_sparse` - Tells `russell` that the local linear solvers are available locally + +### Arch Linux + +#### 1. Default -Then: +Install some dependencies: ```bash -cargo add gemlab + pacman -Syu blas-openblas python-matplotlib suitesparse ``` -### Details +#### 2. Alternative -This crates depends on `russell_lab` and, hence, needs some external libraries. See the [installation of required dependencies](https://github.com/cpmech/russell) on `russell_lab`. +Check [Russell (Rust Scientific Library)](https://github.com/cpmech/russell) out for detailed instructions on how to install the optional dependencies. -### Setting Cargo.toml +### Debian/Ubuntu Linux -[![Crates.io](https://img.shields.io/crates/v/gemlab.svg)](https://crates.io/crates/gemlab) +#### 1. Default -๐Ÿ‘† Check the crate version and update your Cargo.toml accordingly: +Install some dependencies: -```toml -[dependencies] -gemlab = "*" +```bash +sudo apt-get install -y liblapacke-dev libopenblas-dev libsuitesparse-dev python3-matplotlib ``` +#### 2. Alternative + +Check [Russell (Rust Scientific Library)](https://github.com/cpmech/russell) out for detailed instructions on how to install the optional dependencies. + ## Examples diff --git a/data/figures/gauss_points_hex_14.svg b/data/figures/gauss_points_hex_14.svg new file mode 100644 index 0000000..5010d3a --- /dev/null +++ b/data/figures/gauss_points_hex_14.svg @@ -0,0 +1,1248 @@ + + + + + + + + 2026-04-09T15:51:49.083904 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_hex_27.svg b/data/figures/gauss_points_hex_27.svg new file mode 100644 index 0000000..c381df4 --- /dev/null +++ b/data/figures/gauss_points_hex_27.svg @@ -0,0 +1,1443 @@ + + + + + + + + 2026-04-09T15:51:49.774886 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_hex_6.svg b/data/figures/gauss_points_hex_6.svg new file mode 100644 index 0000000..50876a3 --- /dev/null +++ b/data/figures/gauss_points_hex_6.svg @@ -0,0 +1,1058 @@ + + + + + + + + 2026-04-09T15:51:47.724377 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_hex_64.svg b/data/figures/gauss_points_hex_64.svg new file mode 100644 index 0000000..70e98c3 --- /dev/null +++ b/data/figures/gauss_points_hex_64.svg @@ -0,0 +1,1998 @@ + + + + + + + + 2026-04-09T15:51:50.388771 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_hex_8.svg b/data/figures/gauss_points_hex_8.svg new file mode 100644 index 0000000..7055484 --- /dev/null +++ b/data/figures/gauss_points_hex_8.svg @@ -0,0 +1,1127 @@ + + + + + + + + 2026-04-09T15:51:48.427759 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_qua_1.svg b/data/figures/gauss_points_qua_1.svg new file mode 100644 index 0000000..7855d75 --- /dev/null +++ b/data/figures/gauss_points_qua_1.svg @@ -0,0 +1,716 @@ + + + + + + + + 2026-04-09T15:51:47.718891 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_qua_16.svg b/data/figures/gauss_points_qua_16.svg new file mode 100644 index 0000000..495c463 --- /dev/null +++ b/data/figures/gauss_points_qua_16.svg @@ -0,0 +1,1476 @@ + + + + + + + + 2026-04-09T15:51:49.762398 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_qua_4.svg b/data/figures/gauss_points_qua_4.svg new file mode 100644 index 0000000..f90ca30 --- /dev/null +++ b/data/figures/gauss_points_qua_4.svg @@ -0,0 +1,896 @@ + + + + + + + + 2026-04-09T15:51:48.407790 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_qua_9.svg b/data/figures/gauss_points_qua_9.svg new file mode 100644 index 0000000..a33a59e --- /dev/null +++ b/data/figures/gauss_points_qua_9.svg @@ -0,0 +1,1193 @@ + + + + + + + + 2026-04-09T15:51:49.083525 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tet_1.svg b/data/figures/gauss_points_tet_1.svg new file mode 100644 index 0000000..65f6f06 --- /dev/null +++ b/data/figures/gauss_points_tet_1.svg @@ -0,0 +1,731 @@ + + + + + + + + 2026-04-09T15:51:47.715616 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tet_14.svg b/data/figures/gauss_points_tet_14.svg new file mode 100644 index 0000000..87c14e7 --- /dev/null +++ b/data/figures/gauss_points_tet_14.svg @@ -0,0 +1,989 @@ + + + + + + + + 2026-04-09T15:51:50.357796 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tet_15.svg b/data/figures/gauss_points_tet_15.svg new file mode 100644 index 0000000..0b2e54b --- /dev/null +++ b/data/figures/gauss_points_tet_15.svg @@ -0,0 +1,1004 @@ + + + + + + + + 2026-04-09T15:51:50.862571 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tet_24.svg b/data/figures/gauss_points_tet_24.svg new file mode 100644 index 0000000..3b95882 --- /dev/null +++ b/data/figures/gauss_points_tet_24.svg @@ -0,0 +1,1139 @@ + + + + + + + + 2026-04-09T15:51:51.374250 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tet_4.svg b/data/figures/gauss_points_tet_4.svg new file mode 100644 index 0000000..8afaca8 --- /dev/null +++ b/data/figures/gauss_points_tet_4.svg @@ -0,0 +1,773 @@ + + + + + + + + 2026-04-09T15:51:48.407954 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tet_5.svg b/data/figures/gauss_points_tet_5.svg new file mode 100644 index 0000000..389aed1 --- /dev/null +++ b/data/figures/gauss_points_tet_5.svg @@ -0,0 +1,814 @@ + + + + + + + + 2026-04-09T15:51:49.074277 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tet_8.svg b/data/figures/gauss_points_tet_8.svg new file mode 100644 index 0000000..33b1190 --- /dev/null +++ b/data/figures/gauss_points_tet_8.svg @@ -0,0 +1,868 @@ + + + + + + + + 2026-04-09T15:51:49.752181 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tri_1.svg b/data/figures/gauss_points_tri_1.svg new file mode 100644 index 0000000..9390c6e --- /dev/null +++ b/data/figures/gauss_points_tri_1.svg @@ -0,0 +1,674 @@ + + + + + + + + 2026-04-09T15:51:47.706359 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tri_12.svg b/data/figures/gauss_points_tri_12.svg new file mode 100644 index 0000000..79bae5e --- /dev/null +++ b/data/figures/gauss_points_tri_12.svg @@ -0,0 +1,1129 @@ + + + + + + + + 2026-04-09T15:51:50.854093 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tri_16.svg b/data/figures/gauss_points_tri_16.svg new file mode 100644 index 0000000..790cea1 --- /dev/null +++ b/data/figures/gauss_points_tri_16.svg @@ -0,0 +1,1285 @@ + + + + + + + + 2026-04-09T15:51:51.363818 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tri_3.svg b/data/figures/gauss_points_tri_3.svg new file mode 100644 index 0000000..2525bd7 --- /dev/null +++ b/data/figures/gauss_points_tri_3.svg @@ -0,0 +1,752 @@ + + + + + + + + 2026-04-09T15:51:48.400151 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tri_4.svg b/data/figures/gauss_points_tri_4.svg new file mode 100644 index 0000000..104eac5 --- /dev/null +++ b/data/figures/gauss_points_tri_4.svg @@ -0,0 +1,790 @@ + + + + + + + + 2026-04-09T15:51:49.060696 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tri_6.svg b/data/figures/gauss_points_tri_6.svg new file mode 100644 index 0000000..d544ef3 --- /dev/null +++ b/data/figures/gauss_points_tri_6.svg @@ -0,0 +1,898 @@ + + + + + + + + 2026-04-09T15:51:49.751462 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/data/figures/gauss_points_tri_7.svg b/data/figures/gauss_points_tri_7.svg new file mode 100644 index 0000000..bebfcd9 --- /dev/null +++ b/data/figures/gauss_points_tri_7.svg @@ -0,0 +1,936 @@ + + + + + + + + 2026-04-09T15:51:50.347761 + image/svg+xml + + + Matplotlib v3.6.3, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/integ/gauss.rs b/src/integ/gauss.rs index ec64f5e..d6098e0 100644 --- a/src/integ/gauss.rs +++ b/src/integ/gauss.rs @@ -4,6 +4,43 @@ use russell_lab::Matrix; /// Holds a set of integration points ("Gauss") associated with a GeoClass /// +/// # Illustration +/// +/// ## Tri +/// +/// ![Tri](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tri_1.svg) +/// ![Tri](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tri_3.svg) +/// ![Tri](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tri_4.svg) +/// ![Tri](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tri_6.svg) +/// ![Tri](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tri_7.svg) +/// ![Tri](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tri_12.svg) +/// ![Tri](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tri_16.svg) +/// +/// ## Qua +/// +/// ![Qua](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_qua_1.svg) +/// ![Qua](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_qua_4.svg) +/// ![Qua](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_qua_9.svg) +/// ![Qua](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_qua_16.svg) +/// +/// ## Tet +/// +/// ![Tet](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tet_1.svg) +/// ![Tet](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tet_4.svg) +/// ![Tet](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tet_5.svg) +/// ![Tet](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tet_8.svg) +/// ![Tet](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tet_14.svg) +/// ![Tet](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tet_15.svg) +/// ![Tet](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_tet_24.svg) +/// +/// ## Hex +/// +/// ![Hex](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_hex_6.svg) +/// ![Hex](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_hex_8.svg) +/// ![Hex](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_hex_14.svg) +/// ![Hex](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_hex_27.svg) +/// ![Hex](https://raw.githubusercontent.com/cpmech/gemlab/main/data/figures/gauss_points_hex_64.svg) +/// /// # Reference /// /// 1. Durand R and Farias MM (2014) A local extrapolation method for finite elements, @@ -835,4 +872,215 @@ mod tests { assert_eq!(Gauss::new_or_sized(GeoKind::Qua8, None).unwrap().npoint(), 9); assert_eq!(Gauss::new_or_sized(GeoKind::Qua8, Some(1)).unwrap().npoint(), 1); } + + // Do not delete the code below (used to draw the Gauss points) + + /* + + use plotpy::{Canvas, Plot, Text}; + + #[test] + fn draw_gauss_points_tri() { + for ng in [1, 3, 4, 6, 7, 12, 16] { + let mut plot = Plot::new(); + let mut canvas = Canvas::new(); + canvas + .set_face_color("None") + .set_line_width(2.0) + .draw_polyline(&[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]], true); + plot.add(&canvas); + let color = "#000000"; + let size = 0.01; + let gap = 0.015; + canvas.set_face_color(color).set_edge_color(color); + let mut text1 = Text::new(); + let mut text2 = Text::new(); + text1 + .set_color("#eb0707") + .set_fontsize(12.0) + .set_align_horizontal("center") + .set_align_vertical("bottom"); + text2.set_align_horizontal("center").set_align_vertical("top"); + let gauss = Gauss::new_sized(GeoClass::Tri, ng).unwrap(); + for (i, point) in gauss.data.iter().enumerate() { + text1.draw(point[0], point[1] + gap, &format!("{}", 1 + i)); + text2.draw(point[0], point[1] - gap, &format!("({:.3}, {:.3})", point[0], point[1])); + canvas.draw_circle(point[0], point[1], size); + } + plot.add(&canvas) + .add(&text1) + .add(&text2) + .set_range(-0.1, 1.1, -0.1, 1.1) + .set_rotation_ticks_x(90.0) + .set_ticks_x(0.25, -1.0, "") + .set_equal_axes(true) + .grid_and_labels("$\\xi_1$", "$\\xi_2$") + .set_figure_size_points(500.0, 500.0) + .save(&format!("/tmp/gemlab/gauss_points_tri_{}.svg", ng)) + .unwrap(); + } + } + + #[test] + fn draw_gauss_points_qua() { + const COMPARE_HYPLAS: bool = false; + const HYPLAS4: [[f64; 2]; 4] = [ + [-0.577350269189626, -0.577350269189626], + [-0.577350269189626, 0.577350269189626], + [0.577350269189626, -0.577350269189626], + [0.577350269189626, 0.577350269189626], + ]; + const HYPLAS9: [[f64; 2]; 9] = [ + [-0.774596669241483, -0.774596669241483], + [-0.774596669241483, 0.0], + [-0.774596669241483, 0.774596669241483], + [0.0, -0.774596669241483], + [0.0, 0.0], + [0.0, 0.774596669241483], + [0.774596669241483, -0.774596669241483], + [0.774596669241483, 0.0], + [0.774596669241483, 0.774596669241483], + ]; + for ng in [1, 4, 9, 16] { + let mut plot = Plot::new(); + let mut canvas = Canvas::new(); + canvas + .set_face_color("None") + .set_line_width(2.0) + .draw_rectangle(-1.0, -1.0, 2.0, 2.0); + plot.add(&canvas); + let color = "#000000"; + let size = 0.025; + let gap = 0.03; + canvas.set_face_color(color).set_edge_color(color); + let mut text1 = Text::new(); + let mut text2 = Text::new(); + text1 + .set_color("#eb0707") + .set_fontsize(12.0) + .set_align_horizontal("center") + .set_align_vertical("bottom"); + text2.set_align_horizontal("center").set_align_vertical("top"); + let gauss = Gauss::new_sized(GeoClass::Qua, ng).unwrap(); + for (i, point) in gauss.data.iter().enumerate() { + text1.draw(point[0], point[1] + gap, &format!("{}", 1 + i)); + text2.draw(point[0], point[1] - gap, &format!("({:.3}, {:.3})", point[0], point[1])); + canvas.draw_circle(point[0], point[1], size); + } + if COMPARE_HYPLAS && (ng == 4 || ng == 9) { + canvas.set_face_color("None").set_edge_color("green"); + text1.set_align_vertical("top"); + for i in 0..ng { + if ng == 4 { + canvas.draw_circle(HYPLAS4[i][0], HYPLAS4[i][1], 2.0 * size); + text1.draw(HYPLAS4[i][0], HYPLAS4[i][1] - 2.0 * gap, &format!("{}", 1 + i)); + } else { + canvas.draw_circle(HYPLAS9[i][0], HYPLAS9[i][1], 2.0 * size); + text1.draw(HYPLAS9[i][0], HYPLAS9[i][1] - 2.0 * gap, &format!("{}", 1 + i)); + } + } + } + plot.add(&canvas) + .add(&text1) + .add(&text2) + .set_range(-1.1, 1.1, -1.1, 1.1) + .set_rotation_ticks_x(90.0) + .set_ticks_x(0.25, -1.0, "") + .set_equal_axes(true) + .grid_and_labels("$\\xi_1$", "$\\xi_2$") + .set_figure_size_points(500.0, 500.0) + .save(&format!("/tmp/gemlab/gauss_points_qua_{}.svg", ng)) + .unwrap(); + } + } + + #[test] + fn draw_gauss_points_tet() { + for ng in [1, 4, 5, 8, 14, 15, 24] { + let mut plot = Plot::new(); + let mut canvas = Canvas::new(); + canvas.set_face_color("None").set_line_width(2.0).draw_polyline( + &[ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + [1.0, 0.0, 0.0], + ], + false, + ); + canvas.draw_polyline(&[[0.0, 0.0, 1.0], [0.0, 1.0, 0.0]], false); + plot.add(&canvas); + let mut text1 = Text::new(); + text1 + .set_color("#eb0707") + .set_fontsize(12.0) + .set_bbox(true) + .set_bbox_facecolor("white") + .set_bbox_style("square,pad=0.05") + .set_align_horizontal("center") + .set_align_vertical("center"); + let gauss = Gauss::new_sized(GeoClass::Tet, ng).unwrap(); + for (i, point) in gauss.data.iter().enumerate() { + text1.draw_3d(point[0], point[1], point[2], &format!("{}", 1 + i)); + } + plot.add(&canvas) + .add(&text1) + .set_equal_axes(true) + .set_camera(30.0, 15.0) + .set_labels_3d("$\\xi_1$", "$\\xi_2$", "$\\xi_3$") + .set_figure_size_points(500.0, 500.0) + .save(&format!("/tmp/gemlab/gauss_points_tet_{}.svg", ng)) + .unwrap(); + } + } + + #[test] + fn draw_gauss_points_hex() { + for ng in [6, 8, 14, 27, 64] { + let mut plot = Plot::new(); + let mut canvas = Canvas::new(); + canvas.set_face_color("None").set_line_width(2.0).draw_polyline( + &[ + [-1.0, -1.0, -1.0], + [1.0, -1.0, -1.0], + [1.0, 1.0, -1.0], + [-1.0, 1.0, -1.0], + ], + true, + ); + canvas.draw_polyline( + &[[-1.0, -1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, 1.0], [-1.0, 1.0, 1.0]], + true, + ); + canvas.draw_polyline(&[[-1.0, -1.0, -1.0], [-1.0, -1.0, 1.0]], false); + canvas.draw_polyline(&[[1.0, -1.0, -1.0], [1.0, -1.0, 1.0]], false); + canvas.draw_polyline(&[[1.0, 1.0, -1.0], [1.0, 1.0, 1.0]], false); + canvas.draw_polyline(&[[-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0]], false); + plot.add(&canvas); + let mut text1 = Text::new(); + text1 + .set_color("#eb0707") + .set_fontsize(12.0) + .set_bbox(true) + .set_bbox_facecolor("white") + .set_bbox_style("square,pad=0.05") + .set_align_horizontal("center") + .set_align_vertical("center"); + let gauss = Gauss::new_sized(GeoClass::Hex, ng).unwrap(); + for (i, point) in gauss.data.iter().enumerate() { + text1.draw_3d(point[0], point[1], point[2], &format!("{}", 1 + i)); + } + plot.add(&canvas) + .add(&text1) + .set_equal_axes(true) + .set_camera(30.0, 15.0) + .set_labels_3d("$\\xi_1$", "$\\xi_2$", "$\\xi_3$") + .set_figure_size_points(500.0, 500.0) + .save(&format!("/tmp/gemlab/gauss_points_hex_{}.svg", ng)) + .unwrap(); + } + } + */ } diff --git a/src/mesh/block.rs b/src/mesh/block.rs index be690f3..22b49d0 100644 --- a/src/mesh/block.rs +++ b/src/mesh/block.rs @@ -1809,15 +1809,15 @@ mod tests { let mesh = block.subdivide(GeoKind::Qua4).unwrap(); assert_eq!(mesh.points.len(), 9); assert_eq!(mesh.cells.len(), 4); - array_approx_eq(&mesh.points[0].coords, &[0.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[1].coords, &[1.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[2].coords, &[1.0, 1.0], 1e-15); - array_approx_eq(&mesh.points[3].coords, &[0.0, 1.0], 1e-15); - array_approx_eq(&mesh.points[4].coords, &[3.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[5].coords, &[3.0, 1.0], 1e-15); - array_approx_eq(&mesh.points[6].coords, &[1.0, 3.0], 1e-15); - array_approx_eq(&mesh.points[7].coords, &[0.0, 3.0], 1e-15); - array_approx_eq(&mesh.points[8].coords, &[3.0, 3.0], 1e-15); + array_approx_eq(&mesh.points[0].coords, &[0.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[1].coords, &[1.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[2].coords, &[1.0, 1.0], 1e-14); + array_approx_eq(&mesh.points[3].coords, &[0.0, 1.0], 1e-14); + array_approx_eq(&mesh.points[4].coords, &[3.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[5].coords, &[3.0, 1.0], 1e-14); + array_approx_eq(&mesh.points[6].coords, &[1.0, 3.0], 1e-14); + array_approx_eq(&mesh.points[7].coords, &[0.0, 3.0], 1e-14); + array_approx_eq(&mesh.points[8].coords, &[3.0, 3.0], 1e-14); assert_eq!(&mesh.cells[0].points, &[0, 1, 2, 3]); assert_eq!(&mesh.cells[1].points, &[1, 4, 5, 2]); assert_eq!(&mesh.cells[2].points, &[3, 2, 6, 7]); @@ -1934,7 +1934,7 @@ mod tests { assert_eq!(mesh.points.len(), correct.points.len()); assert_eq!(mesh.cells.len(), correct.cells.len()); for point in &correct.points { - array_approx_eq(&point.coords, &correct.points[point.id].coords, 1e-15); + array_approx_eq(&point.coords, &correct.points[point.id].coords, 1e-14); } for cell in &correct.cells { assert_eq!(cell.points, correct.cells[cell.id].points); @@ -1969,7 +1969,7 @@ mod tests { assert_eq!(mesh.points.len(), correct.points.len()); assert_eq!(mesh.cells.len(), correct.cells.len()); for point in &correct.points { - array_approx_eq(&point.coords, &correct.points[point.id].coords, 1e-15); + array_approx_eq(&point.coords, &correct.points[point.id].coords, 1e-14); } for cell in &correct.cells { assert_eq!(cell.points, correct.cells[cell.id].points); @@ -2104,24 +2104,24 @@ mod tests { mesh.check_all().unwrap(); assert_eq!(mesh.points.len(), 27); assert_eq!(mesh.cells.len(), 8); - array_approx_eq(&mesh.points[0].coords, &[0.0, 0.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[1].coords, &[1.0, 0.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[2].coords, &[1.0, 1.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[3].coords, &[0.0, 1.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[4].coords, &[0.0, 0.0, 1.0], 1e-15); - array_approx_eq(&mesh.points[5].coords, &[1.0, 0.0, 1.0], 1e-15); - array_approx_eq(&mesh.points[6].coords, &[1.0, 1.0, 1.0], 1e-15); - array_approx_eq(&mesh.points[7].coords, &[0.0, 1.0, 1.0], 1e-15); - array_approx_eq(&mesh.points[8].coords, &[3.0, 0.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[9].coords, &[3.0, 1.0, 0.0], 1e-15); - array_approx_eq(&mesh.points[10].coords, &[3.0, 0.0, 1.0], 1e-15); + array_approx_eq(&mesh.points[0].coords, &[0.0, 0.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[1].coords, &[1.0, 0.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[2].coords, &[1.0, 1.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[3].coords, &[0.0, 1.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[4].coords, &[0.0, 0.0, 1.0], 1e-14); + array_approx_eq(&mesh.points[5].coords, &[1.0, 0.0, 1.0], 1e-14); + array_approx_eq(&mesh.points[6].coords, &[1.0, 1.0, 1.0], 1e-14); + array_approx_eq(&mesh.points[7].coords, &[0.0, 1.0, 1.0], 1e-14); + array_approx_eq(&mesh.points[8].coords, &[3.0, 0.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[9].coords, &[3.0, 1.0, 0.0], 1e-14); + array_approx_eq(&mesh.points[10].coords, &[3.0, 0.0, 1.0], 1e-14); array_approx_eq(&mesh.points[11].coords, &[3.0, 1.0, 1.0], 1e-14); // << - array_approx_eq(&mesh.points[22].coords, &[3.0, 0.0, 3.0], 1e-15); - array_approx_eq(&mesh.points[23].coords, &[3.0, 1.0, 3.0], 1e-15); - array_approx_eq(&mesh.points[18].coords, &[0.0, 0.0, 3.0], 1e-15); - array_approx_eq(&mesh.points[19].coords, &[1.0, 0.0, 3.0], 1e-15); - array_approx_eq(&mesh.points[20].coords, &[1.0, 1.0, 3.0], 1e-15); - array_approx_eq(&mesh.points[21].coords, &[0.0, 1.0, 3.0], 1e-15); + array_approx_eq(&mesh.points[22].coords, &[3.0, 0.0, 3.0], 1e-14); + array_approx_eq(&mesh.points[23].coords, &[3.0, 1.0, 3.0], 1e-14); + array_approx_eq(&mesh.points[18].coords, &[0.0, 0.0, 3.0], 1e-14); + array_approx_eq(&mesh.points[19].coords, &[1.0, 0.0, 3.0], 1e-14); + array_approx_eq(&mesh.points[20].coords, &[1.0, 1.0, 3.0], 1e-14); + array_approx_eq(&mesh.points[21].coords, &[0.0, 1.0, 3.0], 1e-14); assert_eq!(&mesh.cells[0].points, &[0, 1, 2, 3, 4, 5, 6, 7]); assert_eq!(&mesh.cells[1].points, &[1, 8, 9, 2, 5, 10, 11, 6]); assert_eq!(&mesh.cells[7].points, &[6, 11, 17, 14, 20, 23, 26, 24]); @@ -2225,7 +2225,7 @@ mod tests { } radius = f64::sqrt(radius); if [0, 3, 7].contains(&point.id) { - approx_eq(radius, block.args_ring.rmin, 1e-15); + approx_eq(radius, block.args_ring.rmin, 1e-14); } if [1, 2, 6].contains(&point.id) { approx_eq(radius, (block.args_ring.rmin + block.args_ring.rmax) / 2.0, 1e-17); @@ -2268,7 +2268,7 @@ mod tests { } radius = f64::sqrt(radius); if [0, 11, 7, 3, 35, 32, 29].contains(&point.id) { - approx_eq(radius, block.args_ring.rmin, 1e-15); + approx_eq(radius, block.args_ring.rmin, 1e-14); } if [1, 5, 9, 2, 30, 33, 28].contains(&point.id) { approx_eq(radius, (block.args_ring.rmin + block.args_ring.rmax) / 2.0, 1e-17); @@ -2315,7 +2315,7 @@ mod tests { } radius = f64::sqrt(radius); if [0, 3, 4, 7, 13, 15, 18, 21, 25].contains(&point.id) { - approx_eq(radius, block.args_ring.rmin, 1e-15); + approx_eq(radius, block.args_ring.rmin, 1e-14); } if [1, 2, 5, 6, 12, 14, 19, 20, 24].contains(&point.id) { approx_eq(radius, (block.args_ring.rmin + block.args_ring.rmax) / 2.0, 1e-17); @@ -2324,17 +2324,17 @@ mod tests { approx_eq(radius, block.args_ring.rmax, 1e-17); } if [0, 1, 2, 3, 8, 9, 12, 13, 16].contains(&point.id) { - approx_eq(point.coords[2], block.args_ring.zmin, 1e-15); + approx_eq(point.coords[2], block.args_ring.zmin, 1e-14); } if [4, 5, 6, 7, 10, 11, 14, 15, 17].contains(&point.id) { approx_eq( point.coords[2], (block.args_ring.zmin + block.args_ring.zmax) / 2.0, - 1e-15, + 1e-14, ); } if [18, 19, 20, 21, 22, 23, 24, 25, 26].contains(&point.id) { - approx_eq(point.coords[2], block.args_ring.zmax, 1e-15); + approx_eq(point.coords[2], block.args_ring.zmax, 1e-14); } } if SAVE_FIGURE { @@ -2370,7 +2370,7 @@ mod tests { } radius = f64::sqrt(radius); if [0, 15, 14, 3, 61, 60, 53, 24, 25, 4, 7, 67, 66, 55, 123, 122, 71, 94].contains(&point.id) { - approx_eq(radius, block.args_ring.rmin, 1e-15); + approx_eq(radius, block.args_ring.rmin, 1e-14); } if [90, 91, 118, 62, 56, 63, 57, 119, 5, 27, 1].contains(&point.id) { approx_eq(radius, (block.args_ring.rmin + block.args_ring.rmax) / 2.0, 1e-17); @@ -2423,7 +2423,7 @@ mod tests { let mesh = block.subdivide(GeoKind::Qua4).unwrap(); for p in [0, 3, 7] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, 0.0]).unwrap(); - approx_eq(d, 1.0, 1e-15); + approx_eq(d, 1.0, 1e-14); } if SAVE_FIGURE { draw( @@ -2441,7 +2441,7 @@ mod tests { let mesh = block.subdivide(GeoKind::Qua4).unwrap(); for p in [0, 3, 7] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, 0.0]).unwrap(); - approx_eq(d, 0.5, 1e-15); + approx_eq(d, 0.5, 1e-14); } if SAVE_FIGURE { draw( @@ -2475,13 +2475,13 @@ mod tests { // side 3 for p in [0, 3, 7, 14, 17] { let d = point_point_distance(&mesh.points[p].coords, &[xc, xc]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // middle nodes for (a, mid, b) in [(3, 6, 2)] { let xmid = (mesh.points[a].coords[0] + mesh.points[b].coords[0]) / 2.0; let ymid = (mesh.points[a].coords[1] + mesh.points[b].coords[1]) / 2.0; - array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-15); + array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-14); } if SAVE_FIGURE { draw( @@ -2526,22 +2526,22 @@ mod tests { // side 0 for p in [0, 4, 1, 10, 8] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, cen_minus]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // side 1 for p in [8, 11, 9, 19, 18] { let d = point_point_distance(&mesh.points[p].coords, &[cen_plus, 0.0]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // side 2 for p in [14, 16, 13, 20, 18] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, cen_plus]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // side 3 for p in [0, 7, 3, 17, 14] { let d = point_point_distance(&mesh.points[p].coords, &[cen_minus, 0.0]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // center vertical line for p in [5, 2, 15] { @@ -2555,7 +2555,7 @@ mod tests { for (a, mid, b) in [(3, 6, 2), (2, 12, 9), (1, 5, 2), (2, 15, 13)] { let xmid = (mesh.points[a].coords[0] + mesh.points[b].coords[0]) / 2.0; let ymid = (mesh.points[a].coords[1] + mesh.points[b].coords[1]) / 2.0; - array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-15); + array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-14); } if SAVE_FIGURE { draw( @@ -2617,7 +2617,7 @@ mod tests { ] { let x = &[mesh.points[p].coords[0], mesh.points[p].coords[1]]; let d = point_point_distance(x, &[cen_minus, 0.0]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // face 1 for p in [ @@ -2629,7 +2629,7 @@ mod tests { ] { let x = &[mesh.points[p].coords[0], mesh.points[p].coords[1]]; let d = point_point_distance(x, &[cen_plus, 0.0]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // face 2 for p in [ @@ -2641,7 +2641,7 @@ mod tests { ] { let x = &[mesh.points[p].coords[0], mesh.points[p].coords[1]]; let d = point_point_distance(x, &[0.0, cen_minus]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // face 3 for p in [ @@ -2653,7 +2653,7 @@ mod tests { ] { let x = &[mesh.points[p].coords[0], mesh.points[p].coords[1]]; let d = point_point_distance(x, &[0.0, cen_plus]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // vertical line at center for p in [2, 18, 6, 61, 53] { @@ -2678,7 +2678,7 @@ mod tests { let xmid = (mesh.points[a].coords[0] + mesh.points[b].coords[0]) / 2.0; let ymid = (mesh.points[a].coords[1] + mesh.points[b].coords[1]) / 2.0; let zmid = (mesh.points[a].coords[2] + mesh.points[b].coords[2]) / 2.0; - array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid, zmid], 1e-15); + array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid, zmid], 1e-14); } if SAVE_FIGURE { draw( @@ -2751,7 +2751,7 @@ mod tests { ] { let x = &[mesh.points[p].coords[1], mesh.points[p].coords[2]]; let d = point_point_distance(x, &[0.0, cen_minus]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // face 3 for p in [ @@ -2763,7 +2763,7 @@ mod tests { ] { let x = &[mesh.points[p].coords[0], mesh.points[p].coords[2]]; let d = point_point_distance(x, &[0.0, cen_plus]).unwrap(); - approx_eq(d, r, 1e-15); + approx_eq(d, r, 1e-14); } // vertical line at center for p in [2, 18, 6, 61, 53] { @@ -2782,7 +2782,7 @@ mod tests { let xmid = (mesh.points[a].coords[0] + mesh.points[b].coords[0]) / 2.0; let ymid = (mesh.points[a].coords[1] + mesh.points[b].coords[1]) / 2.0; let zmid = (mesh.points[a].coords[2] + mesh.points[b].coords[2]) / 2.0; - array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid, zmid], 1e-15); + array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid, zmid], 1e-14); } if SAVE_FIGURE { draw( @@ -2816,11 +2816,11 @@ mod tests { let mesh = block.subdivide(GeoKind::Qua4).unwrap(); for p in [6, 7, 11, 15] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, 0.0]).unwrap(); - approx_eq(d, 6.0, 1e-15); + approx_eq(d, 6.0, 1e-14); } for p in [13, 12, 14] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, 0.0]).unwrap(); - approx_eq(d, 6.0, 1e-15); + approx_eq(d, 6.0, 1e-14); } if SAVE_FIGURE { draw( @@ -2846,12 +2846,12 @@ mod tests { let mesh = block.subdivide(GeoKind::Qua8).unwrap(); for p in [13, 16, 14, 27, 26, 38, 37] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, 0.0]).unwrap(); - approx_eq(d, radius, 1e-15); + approx_eq(d, radius, 1e-14); } for (a, mid, b) in [(23, 28, 26), (9, 17, 14)] { let xmid = (mesh.points[a].coords[0] + mesh.points[b].coords[0]) / 2.0; let ymid = (mesh.points[a].coords[1] + mesh.points[b].coords[1]) / 2.0; - array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-15); + array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-14); } if SAVE_FIGURE { draw( @@ -2877,12 +2877,12 @@ mod tests { let mesh = block.subdivide(GeoKind::Qua9).unwrap(); for p in [15, 18, 16, 32, 31, 46, 45] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, 0.0]).unwrap(); - approx_eq(d, radius, 1e-15); + approx_eq(d, radius, 1e-14); } for (a, mid, b) in [(42, 48, 46), (27, 33, 31), (28, 34, 32), (10, 19, 16), (12, 20, 18)] { let xmid = (mesh.points[a].coords[0] + mesh.points[b].coords[0]) / 2.0; let ymid = (mesh.points[a].coords[1] + mesh.points[b].coords[1]) / 2.0; - array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-15); + array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-14); } if SAVE_FIGURE { draw( @@ -2923,7 +2923,7 @@ mod tests { let xd = mesh.points[a].coords[0] + 2.0 * lx / 3.0; let yd = mesh.points[a].coords[1] + 2.0 * ly / 3.0; array_approx_eq(&mesh.points[c].coords, &[xc, yc], 1e-14); - array_approx_eq(&mesh.points[d].coords, &[xd, yd], 1e-15); + array_approx_eq(&mesh.points[d].coords, &[xd, yd], 1e-14); } if SAVE_FIGURE { draw( @@ -2950,7 +2950,7 @@ mod tests { for (a, mid, b) in [(82, 92, 90), (54, 64, 62), (20, 34, 32)] { let xmid = (mesh.points[a].coords[0] + mesh.points[b].coords[0]) / 2.0; let ymid = (mesh.points[a].coords[1] + mesh.points[b].coords[1]) / 2.0; - array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-15); + array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid], 1e-14); } if SAVE_FIGURE { draw( @@ -2990,13 +2990,13 @@ mod tests { 22, 28, 23, 48, 45, // z-max ] { let d = point_point_distance(&mesh.points[p].coords, &[0.0, 0.0, mesh.points[p].coords[2]]).unwrap(); - approx_eq(d, radius, 1e-15); + approx_eq(d, radius, 1e-14); } for (a, mid, b) in [(2, 26, 21), (6, 29, 23)] { let xmid = (mesh.points[a].coords[0] + mesh.points[b].coords[0]) / 2.0; let ymid = (mesh.points[a].coords[1] + mesh.points[b].coords[1]) / 2.0; let zmid = (mesh.points[a].coords[2] + mesh.points[b].coords[2]) / 2.0; - array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid, zmid], 1e-15); + array_approx_eq(&mesh.points[mid].coords, &[xmid, ymid, zmid], 1e-14); } if SAVE_FIGURE { draw( diff --git a/src/mesh/edges.rs b/src/mesh/edges.rs index 091c200..abce706 100644 --- a/src/mesh/edges.rs +++ b/src/mesh/edges.rs @@ -243,6 +243,38 @@ impl<'a> Edges<'a> { (path, points) } + + /// Returns a sorted list of all unique points in the collection of edges + /// + /// # Returns + /// + /// A vector containing all unique point IDs from all edges in `self.all`, sorted in ascending order. + /// + /// # Examples + /// + /// ``` + /// use gemlab::mesh::{Edge, Edges}; + /// use gemlab::shapes::GeoKind; + /// + /// let e1 = Edge { kind: GeoKind::Lin2, points: vec![1, 2], marker: 0 }; + /// let e2 = Edge { kind: GeoKind::Lin2, points: vec![2, 3], marker: 0 }; + /// let e3 = Edge { kind: GeoKind::Lin3, points: vec![3, 4, 5], marker: 0 }; + /// let edges = Edges { all: vec![&e1, &e2, &e3] }; + /// + /// let points = edges.all_points(); + /// assert_eq!(points, vec![1, 2, 3, 4, 5]); + /// ``` + pub fn all_points(&self) -> Vec { + let mut points_set = HashSet::new(); + for edge in &self.all { + for &point in &edge.points { + points_set.insert(point); + } + } + let mut points: Vec<_> = points_set.into_iter().collect(); + points.sort(); + points + } } impl fmt::Display for Edge { @@ -586,4 +618,41 @@ mod tests { println!("{:?}", edges.any_path()); assert_eq!(edges.any_path(), (vec![0, 1], vec![3, 14, 6, 13, 2, 28, 21, 27, 18])); } + + #[test] + fn all_points_works() { + // Test with Lin2 edges + let all = generate_sample_lin2(); + + // Empty list of edges + let empty = Edges { all: vec![] }; + assert!(empty.all_points().is_empty()); + + // Single edge + let single = Edges { all: vec![&all[9]] }; + assert_eq!(single.all_points(), vec![3, 4]); + + // Multiple edges with some shared points + let multiple = Edges { + all: vec![&all[0], &all[1], &all[2]], + }; + assert_eq!(multiple.all_points(), vec![2, 3, 5, 10]); + + // All edges (should get all unique points, sorted) + let all_edges = Edges { + all: vec![ + &all[0], &all[1], &all[2], &all[3], &all[4], &all[5], &all[6], &all[7], &all[8], &all[9], &all[10], + &all[11], &all[12], &all[13], + ], + }; + assert_eq!(all_edges.all_points(), vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 10]); + + // Test with Lin3 edges + let all_lin3 = generate_sample_lin3(); + let edges_lin3 = Edges { + all: vec![&all_lin3[0], &all_lin3[1]], + }; + // Points include middle nodes (100, 101) + assert_eq!(edges_lin3.all_points(), vec![2, 3, 5, 100, 101]); + } } diff --git a/src/mesh/faces.rs b/src/mesh/faces.rs index 4a6aff8..cd8ccc1 100644 --- a/src/mesh/faces.rs +++ b/src/mesh/faces.rs @@ -59,6 +59,39 @@ pub type MapFaceToCells = HashMap>; /// Relates a point id to a unique set of FaceKey pub type MapPointToFaces = HashMap>; +impl<'a> Faces<'a> { + /// Returns a sorted list of all unique points in the collection of faces + /// + /// # Returns + /// + /// A vector containing all unique point IDs from all faces in `self.all`, sorted in ascending order. + /// + /// # Examples + /// + /// ``` + /// use gemlab::mesh::{Face, Faces}; + /// use gemlab::shapes::GeoKind; + /// + /// let f1 = Face { kind: GeoKind::Tri3, points: vec![1, 2, 3], marker: 0 }; + /// let f2 = Face { kind: GeoKind::Tri3, points: vec![2, 3, 4], marker: 0 }; + /// let faces = Faces { all: vec![&f1, &f2] }; + /// + /// let points = faces.all_points(); + /// assert_eq!(points, vec![1, 2, 3, 4]); + /// ``` + pub fn all_points(&self) -> Vec { + let mut points_set = HashSet::new(); + for face in &self.all { + for &point in &face.points { + points_set.insert(point); + } + } + let mut points: Vec<_> = points_set.into_iter().collect(); + points.sort(); + points + } +} + impl AsCell for Face { fn kind(&self) -> GeoKind { self.kind @@ -96,3 +129,63 @@ impl<'a> fmt::Display for Faces<'a> { Ok(()) } } + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#[cfg(test)] +mod tests { + use super::{Face, Faces}; + use crate::shapes::GeoKind; + + #[rustfmt::skip] + fn generate_sample_1() -> Vec { + vec![ + Face { kind: GeoKind::Tri3, points: vec![1, 2, 3], marker: 0 }, + Face { kind: GeoKind::Tri3, points: vec![4, 5, 6], marker: 0 }, + ] + } + + #[rustfmt::skip] + fn generate_sample_2() -> Vec { + vec![ + Face { kind: GeoKind::Qua4, points: vec![1, 2, 3, 4], marker: 0 }, + Face { kind: GeoKind::Qua4, points: vec![2, 5, 6, 3], marker: 0 }, + Face { kind: GeoKind::Tri6, points: vec![7, 8, 9, 10, 11, 12], marker: 0 }, + ] + } + + #[test] + fn all_points_works() { + // Empty list of faces + let empty = Faces { all: vec![] }; + assert!(empty.all_points().is_empty()); + + // Single face (Tri3) + let all = generate_sample_1(); + let single = Faces { all: vec![&all[0]] }; + assert_eq!(single.all_points(), vec![1, 2, 3]); + + // Two faces with no shared points + let two_faces = Faces { + all: vec![&all[0], &all[1]], + }; + assert_eq!(two_faces.all_points(), vec![1, 2, 3, 4, 5, 6]); + + // Mixed face types with shared points + let all2 = generate_sample_2(); + let mixed = Faces { + all: vec![&all2[0], &all2[1]], + }; + assert_eq!(mixed.all_points(), vec![1, 2, 3, 4, 5, 6]); + + // Higher-order face (Tri6 with middle nodes) + let higher_order = Faces { all: vec![&all2[2]] }; + assert_eq!(higher_order.all_points(), vec![7, 8, 9, 10, 11, 12]); + + // All faces together + let all_faces = Faces { + all: vec![&all2[0], &all2[1], &all2[2]], + }; + assert_eq!(all_faces.all_points(), vec![1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]); + } +} diff --git a/src/mesh/features.rs b/src/mesh/features.rs index 6cb859f..7625020 100644 --- a/src/mesh/features.rs +++ b/src/mesh/features.rs @@ -348,6 +348,36 @@ impl<'a> Features<'a> { ids } + /// Returns cells sharing a given (2D) edge along with the edge points + /// + /// Returns a **sorted** (by Cell ID) list of tuples containing: + /// * Cell ID + /// * Vector of point IDs on the edge + /// + /// # Examples + /// + /// ``` + /// use gemlab::mesh::{Features, Samples}; + /// + /// // 1.0 4-----------3 + /// // | | + /// // | [1] | + /// // | | + /// // 0.0 0-----------1-----------2-----------5 + /// // [0] [2] + /// let mesh = Samples::mixed_shapes_2d(); + /// let feat = Features::new(&mesh, false); + /// let edge = feat.get_edge(1, 2); + /// let pairs = feat.get_cells_pairs_via_2d_edge(&edge); + /// assert_eq!(pairs, vec![(1, vec![2, 1])]); + /// ``` + pub fn get_cells_pairs_via_2d_edge(&self, edge: &Edge) -> Vec<(CellId, Vec)> { + let cells = self.all_2d_edges.get(&edge.key()).expect("cannot find 2D edge"); + let mut pairs: Vec<_> = cells.iter().map(|c| (c.0, edge.points.clone())).collect(); + pairs.sort_by_key(|p| p.0); + pairs + } + /// Returns cells sharing a given face /// /// Returns a **sorted** list of Cell IDs @@ -358,24 +388,118 @@ impl<'a> Features<'a> { ids } + /// Returns cells sharing a given face along with the face points + /// + /// Returns a **sorted** (by Cell ID) list of tuples containing: + /// * Cell ID + /// * Vector of point IDs on the face + /// + /// # Examples + /// + /// ``` + /// use gemlab::mesh::{Features, Samples}; + /// + /// let mesh = Samples::one_tet4(); + /// let feat = Features::new(&mesh, false); + /// let face = feat.get_face(0, 1, 3, usize::MAX); + /// let pairs = feat.get_cells_pairs_via_face(&face); + /// assert_eq!(pairs, vec![(0, vec![0, 1, 3])]); + /// ``` + pub fn get_cells_pairs_via_face(&self, face: &Face) -> Vec<(CellId, Vec)> { + let cells = self.all_faces.get(&face.key()).expect("cannot find face"); + let mut pairs: Vec<_> = cells.iter().map(|c| (c.0, face.points.clone())).collect(); + pairs.sort_by_key(|p| p.0); + pairs + } + /// Returns many cells sharing a given (2D) edge /// /// Returns a **sorted** list of Cell IDs + /// + /// **Note:** If a cell shares multiple edges from the collection, it will appear multiple times. pub fn get_cells_via_2d_edges(&self, edges: &Edges) -> Vec { let mut ids: Vec<_> = edges.all.iter().flat_map(|e| self.get_cells_via_2d_edge(e)).collect(); ids.sort(); ids } + /// Returns many cells sharing given (2D) edges along with the edge points + /// + /// Returns a **sorted** (by Cell ID) list of tuples containing: + /// * Cell ID + /// * Vector of point IDs on the edge + /// + /// **Note:** If a cell shares multiple edges from the collection, it will appear multiple times. + /// + /// # Examples + /// + /// ``` + /// use gemlab::mesh::{Edges, Features, Samples}; + /// + /// // 1.0 4-----------3 + /// // | | + /// // | [1] | + /// // | | + /// // 0.0 0-----------1-----------2-----------5 + /// // [0] [2] + /// let mesh = Samples::mixed_shapes_2d(); + /// let feat = Features::new(&mesh, false); + /// let edge = feat.get_edge(1, 2); + /// let edges = Edges { all: vec![&edge] }; + /// let pairs = feat.get_cells_pairs_via_2d_edges(&edges); + /// assert_eq!(pairs, vec![(1, vec![2, 1])]); + /// ``` + pub fn get_cells_pairs_via_2d_edges(&self, edges: &Edges) -> Vec<(CellId, Vec)> { + let mut pairs: Vec<_> = edges + .all + .iter() + .flat_map(|e| self.get_cells_pairs_via_2d_edge(e)) + .collect(); + pairs.sort_by_key(|p| p.0); + pairs + } + /// Returns many cells sharing a given face /// /// Returns a **sorted** list of Cell IDs + /// + /// **Note:** If a cell shares multiple faces from the collection, it will appear multiple times. pub fn get_cells_via_faces(&self, faces: &Faces) -> Vec { let mut ids: Vec<_> = faces.all.iter().flat_map(|f| self.get_cells_via_face(f)).collect(); ids.sort(); ids } + /// Returns many cells sharing given faces along with the face points + /// + /// Returns a **sorted** (by Cell ID) list of tuples containing: + /// * Cell ID + /// * Vector of point IDs on the face + /// + /// **Note:** If a cell shares multiple faces from the collection, it will appear multiple times. + /// + /// # Examples + /// + /// ``` + /// use gemlab::mesh::{Faces, Features, Samples}; + /// + /// let mesh = Samples::one_tet4(); + /// let feat = Features::new(&mesh, false); + /// let face = feat.get_face(0, 1, 3, usize::MAX); + /// let faces = Faces { all: vec![&face] }; + /// let pairs = feat.get_cells_pairs_via_faces(&faces); + /// assert_eq!(pairs, vec![(0, vec![0, 1, 3])]); + /// ``` + pub fn get_cells_pairs_via_faces(&self, faces: &Faces) -> Vec<(CellId, Vec)> { + let mut pairs: Vec<_> = faces + .all + .iter() + .flat_map(|f| self.get_cells_pairs_via_face(f)) + .collect(); + pairs.sort_by_key(|p| p.0); + pairs + } + /// Returns all points (sorted) on a set of (2D) edges pub fn get_points_via_2d_edges(&self, edges: &Edges) -> Vec { let mut point_ids: Vec<_> = edges.all.iter().flat_map(|e| e.points.clone()).collect(); @@ -544,7 +668,7 @@ impl<'a> Features<'a> { /// Searches edges with a given marker /// /// Returns edges sorted by their edge keys - pub fn search_marked_edges(&self, marker: i32) -> Edges { + pub fn search_marked_edges(&self, marker: i32) -> Edges<'_> { let mut all: Vec<_> = self.edges.values().filter(|edge| edge.marker == marker).collect(); all.sort_by_key(|edge| edge.key()); Edges { all } @@ -553,7 +677,7 @@ impl<'a> Features<'a> { /// Searches faces with a given marker /// /// Returns faces sorted by their face keys - pub fn search_marked_faces(&self, marker: i32) -> Faces { + pub fn search_marked_faces(&self, marker: i32) -> Faces<'_> { let mut all: Vec<_> = self.faces.values().filter(|face| face.marker == marker).collect(); all.sort_by_key(|face| face.key()); Faces { all } @@ -783,7 +907,7 @@ impl<'a> Features<'a> { /// /// * If at least one point has been found, returns an array such that the edge keys are **sorted** /// * Otherwise, returns an error - pub fn search_edges(&self, at: At, filter: F) -> Result + pub fn search_edges(&self, at: At, filter: F) -> Result, StrError> where F: FnMut(&[f64]) -> bool, { @@ -816,7 +940,7 @@ impl<'a> Features<'a> { /// /// * If at least one point has been found, returns an array such that the face keys are **sorted** /// * Otherwise, returns an error - pub fn search_faces(&self, at: At, filter: F) -> Result + pub fn search_faces(&self, at: At, filter: F) -> Result, StrError> where F: FnMut(&[f64]) -> bool, { @@ -850,7 +974,7 @@ impl<'a> Features<'a> { /// * Returns edges sorted by keys /// * **Warning** Every `At` in the `ats` must generate at least one edge, /// otherwise an error will occur. - pub fn search_many_edges(&self, ats: &[At], mut filter: F) -> Result + pub fn search_many_edges(&self, ats: &[At], mut filter: F) -> Result, StrError> where F: FnMut(&[f64]) -> bool, { @@ -892,7 +1016,7 @@ impl<'a> Features<'a> { /// * Returns faces sorted by keys /// * **Warning** Every `At` in the `ats` must generate at least one face, /// otherwise an error will occur. - pub fn search_many_faces(&self, ats: &[At], mut filter: F) -> Result + pub fn search_many_faces(&self, ats: &[At], mut filter: F) -> Result, StrError> where F: FnMut(&[f64]) -> bool, { @@ -918,6 +1042,88 @@ impl<'a> Features<'a> { Err(e) => Err(e), } } + + /// Search points, cells, and edges in 2D + /// + /// Returns `(point_ids, cell_ids, edges)` + /// + /// # Arguments + /// + /// * `cells_by_points` -- If `true`, all cells connected to the found points will be returned, + /// even if they are not connected to the found edges. If `false`, only cells connected to the found + /// edges will be returned, i.e., `cells_by_points = false` corresponds to `cells_by_edges = true`. + /// * `at` -- the constraint to find **points** + /// * `filter` -- function `fn(x) -> bool` to impose further constraints on the point coordinates + /// + /// Note: Unless `extract_all` was set to `true` when allocating this instance, only boundary edges + /// will be found, affecting the resulting cells as well. + pub fn search_features_2d( + &self, + cells_by_points: bool, + at: At, + mut filter: F, + ) -> Result<(Vec, Vec, Edges<'_>), StrError> + where + F: FnMut(&[f64]) -> bool, + { + // Check + if self.mesh.ndim != 2 { + return Err("search_features_2d works in 2D only"); + } + + // Search all points constrained by "at" and "filter" + let point_ids = self.search_point_ids(at, &mut filter)?; + + // Allocate structures to hold the results + let mut edge_keys_set = HashSet::new(); + let mut cell_ids_set = HashSet::new(); + + // Loop over point ids + for point_id in &point_ids { + // Loop over edges connected to points (using point_to_edges) + if let Some(edge_keys) = self.point_to_edges.get(&point_id) { + for edge_key in edge_keys { + // Accept edge when at least two edge points validate "At" + let accept_edge = if point_ids.contains(&edge_key.0) && point_ids.contains(&edge_key.1) { + edge_keys_set.insert(*edge_key); + true + } else { + false + }; + if cells_by_points || accept_edge { + // Loop over cells connected to edges (using all_2d_edges) + if let Some(cells) = self.all_2d_edges.get(edge_key) { + for &(cell_id, _local_edge_index) in cells { + // Insert a CellId into the HashSet + cell_ids_set.insert(cell_id); + } + } + } + } + } + } + + // Sort edge keys and prepare the Edges structure + let mut edge_keys: Vec<_> = edge_keys_set.into_iter().collect(); + edge_keys.sort(); + let mut edges = Edges { + all: Vec::with_capacity(edge_keys.len()), + }; + for edge_key in edge_keys { + let edge = self + .edges + .get(&edge_key) + .ok_or("INTERNAL ERROR: features.edges data is inconsistent")?; + edges.all.push(edge); + } + + // Sort cell ids + let mut cell_ids: Vec<_> = cell_ids_set.into_iter().collect(); + cell_ids.sort(); + + // Return the results + Ok((point_ids, cell_ids, edges)) + } } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -925,7 +1131,7 @@ impl<'a> Features<'a> { #[cfg(test)] mod tests { use super::{Edge, Edges, Face, Faces, Features}; - use crate::mesh::{At, Draw, Samples}; + use crate::mesh::{At, Draw, Mesh, Samples}; use crate::shapes::GeoKind; use crate::util::any_x; use plotpy::Plot; @@ -988,6 +1194,7 @@ mod tests { assert_eq!(edge.key(), (1, 2)); assert_eq!(face.key(), (0, 1, 3, usize::MAX)); assert_eq!(feat.get_cells_via_face(&face), &[0]); + assert_eq!(feat.get_cells_pairs_via_face(&face), vec![(0, vec![0, 1, 3])]); } #[test] @@ -2163,6 +2370,8 @@ mod tests { assert_eq!(feat.get_cells_via_2d_edge(&edge), &[1]); assert_eq!(feat.get_cells_via_2d_edges(&edges), &[1]); assert_eq!(feat.get_points_via_2d_edges(&edges), &[1, 2]); + assert_eq!(feat.get_cells_pairs_via_2d_edge(&edge), vec![(1, vec![2, 1])]); + assert_eq!(feat.get_cells_pairs_via_2d_edges(&edges), vec![(1, vec![2, 1])]); } #[test] @@ -2191,6 +2400,18 @@ mod tests { assert_eq!(feat.get_cells_via_2d_edge(&edge_b), &[1, 3]); assert_eq!(feat.get_cells_via_2d_edges(&edges), &[0, 1, 2, 3]); assert_eq!(feat.get_points_via_2d_edges(&edges), &[2, 3, 5]); + assert_eq!( + feat.get_cells_pairs_via_2d_edge(&edge_a), + vec![(0, vec![3, 2]), (2, vec![3, 2])] + ); + assert_eq!( + feat.get_cells_pairs_via_2d_edge(&edge_b), + vec![(1, vec![2, 5]), (3, vec![2, 5])] + ); + assert_eq!( + feat.get_cells_pairs_via_2d_edges(&edges), + vec![(0, vec![3, 2]), (1, vec![2, 5]), (2, vec![3, 2]), (3, vec![2, 5])] + ); } #[test] @@ -2219,6 +2440,12 @@ mod tests { assert_eq!(feat.get_cells_via_face(&face_b), &[1]); assert_eq!(feat.get_cells_via_faces(&faces), &[0, 1]); assert_eq!(feat.get_points_via_faces(&faces), &[2, 3, 6, 7]); + assert_eq!(feat.get_cells_pairs_via_face(&face_a), vec![(0, vec![2, 3, 7, 6])]); + assert_eq!(feat.get_cells_pairs_via_face(&face_b), vec![(1, vec![2, 6, 3])]); + assert_eq!( + feat.get_cells_pairs_via_faces(&faces), + vec![(0, vec![2, 3, 7, 6]), (1, vec![2, 6, 3])] + ); } #[test] @@ -2234,5 +2461,278 @@ mod tests { assert_eq!(feat.get_cells_via_face(&face_b), &[4, 6]); assert_eq!(feat.get_cells_via_faces(&faces), &[0, 2, 4, 6]); assert_eq!(feat.get_points_via_faces(&faces), &[2, 3, 6, 7, 20, 21]); + assert_eq!( + feat.get_cells_pairs_via_face(&face_a), + vec![(0, vec![2, 3, 7, 6]), (2, vec![2, 3, 7, 6])] + ); + assert_eq!( + feat.get_cells_pairs_via_face(&face_b), + vec![(4, vec![6, 7, 21, 20]), (6, vec![6, 7, 21, 20])] + ); + assert_eq!( + feat.get_cells_pairs_via_faces(&faces), + vec![ + (0, vec![2, 3, 7, 6]), + (2, vec![2, 3, 7, 6]), + (4, vec![6, 7, 21, 20]), + (6, vec![6, 7, 21, 20]) + ] + ); + } + + #[test] + fn search_features_2d_works_1() { + let mesh = Samples::block_2d_four_qua4(); + + let feat_bry = Features::new(&mesh, false); // no internal edges + let feat_int = Features::new(&mesh, true); // all edges including internal edges + + // cells_by_points ----------------------------------------------------------------------- + let cells_by_points = true; + + // points along x=1.0 (but the internal is not included) + let (point_ids, cell_ids, edges) = feat_bry + .search_features_2d(cells_by_points, At::X(1.0), |_| true) + .unwrap(); + assert_eq!(point_ids, vec![1, 6]); + assert_eq!(cell_ids, vec![0, 1, 2, 3]); + assert_eq!(edges.all.len(), 0); + + // points along x=1.0 with y<=1.0 (but the internal is not included) + let (point_ids, cell_ids, edges) = feat_bry + .search_features_2d(cells_by_points, At::X(1.0), |x| x[1] <= 1.0) + .unwrap(); + assert_eq!(point_ids, vec![1]); + assert_eq!(cell_ids, vec![0, 1]); + assert_eq!(edges.all.len(), 0); + + // points along x=1.0 (with internal edges) + let (point_ids, cell_ids, edges) = feat_int + .search_features_2d(cells_by_points, At::X(1.0), |_| true) + .unwrap(); + assert_eq!(point_ids, vec![1, 2, 6]); + assert_eq!(cell_ids, vec![0, 1, 2, 3]); + assert_eq!( + edges.all.iter().map(|e| e.key()).collect::>(), + vec![(1, 2), (2, 6)] + ); + + // points along x=1.0 with y<=1.0 (with internal edges) + let (point_ids, cell_ids, edges) = feat_int + .search_features_2d(cells_by_points, At::X(1.0), |x| x[1] <= 1.0) + .unwrap(); + assert_eq!(point_ids, vec![1, 2]); + assert_eq!(cell_ids, vec![0, 1, 2, 3]); // cells # 2 and 3 are selected because they share point 2 + assert_eq!(edges.all.iter().map(|e| e.key()).collect::>(), vec![(1, 2)]); + + // cells_by_edges ----------------------------------------------------------------------- + let cells_by_points = false; // only cells attached to edges are returned + + // points along x=1.0 (but the internal is not included) + let (point_ids, cell_ids, edges) = feat_bry + .search_features_2d(cells_by_points, At::X(1.0), |_| true) + .unwrap(); + assert_eq!(point_ids, vec![1, 6]); + assert_eq!(cell_ids.len(), 0); // no cells attached to internal edges are found (because internal edges are not included) + assert_eq!(edges.all.len(), 0); + + // points along x=1.0 with y<=1.0 (but the internal is not included) + let (point_ids, cell_ids, edges) = feat_bry + .search_features_2d(cells_by_points, At::X(1.0), |x| x[1] <= 1.0) + .unwrap(); + assert_eq!(point_ids, vec![1]); + assert_eq!(cell_ids.len(), 0); // no cells attached to internal edges are found (because internal edges are not included) + assert_eq!(edges.all.len(), 0); + + // points along x=1.0 (with internal edges) + let (point_ids, cell_ids, edges) = feat_int + .search_features_2d(cells_by_points, At::X(1.0), |_| true) + .unwrap(); + assert_eq!(point_ids, vec![1, 2, 6]); + assert_eq!(cell_ids, vec![0, 1, 2, 3]); + assert_eq!( + edges.all.iter().map(|e| e.key()).collect::>(), + vec![(1, 2), (2, 6)] + ); + + // points along x=1.0 with y<=1.0 (with internal edges) + let (point_ids, cell_ids, edges) = feat_int + .search_features_2d(cells_by_points, At::X(1.0), |x| x[1] <= 1.0) + .unwrap(); + assert_eq!(point_ids, vec![1, 2]); + assert_eq!(cell_ids, vec![0, 1]); // cells # 2 and 3 are not included because they do not share edge (1, 2) + assert_eq!(edges.all.iter().map(|e| e.key()).collect::>(), vec![(1, 2)]); + } + + #[test] + fn search_features_2d_works_2() { + // L-shaped domain discretized with unstructured triangular elements + let mesh = Mesh::from_text( + r" +# header +# ndim npoint ncell nmarked_edge nmarked_face +2 18 22 12 0 + +# points +# id marker x y {z} +0 1 0.0 0.0 +1 1 2.0 0.0 +2 1 2.0 1.0 +3 1 1.0 1.0 +4 4 1.0 2.0 +5 4 0.0 2.0 +6 1 1.0 0.0 +7 1 0.0 1.0 +8 1 0.5 0.0 +9 0 1.0 0.5 +10 0 1.5 0.25 +11 1 2.0 0.5 +12 0 0.25 0.5 +13 0 0.5 0.9375 +14 0 0.68359375 1.5 +15 4 0.5 2.0 +16 0 1.5 0.7916666666666667 +17 1 0.0 1.5 + +# cells +# id marker kind points +0 1 tri3 0 12 7 +1 1 tri3 6 1 10 +2 1 tri3 8 6 9 +3 1 tri3 13 14 17 +4 1 tri3 16 9 10 +5 1 tri3 10 11 16 +6 1 tri3 6 10 9 +7 1 tri3 8 9 12 +8 1 tri3 12 9 13 +9 1 tri3 14 13 3 +10 1 tri3 2 16 11 +11 1 tri3 13 9 3 +12 1 tri3 0 8 12 +13 1 tri3 10 1 11 +14 1 tri3 3 4 14 +15 1 tri3 7 12 13 +16 1 tri3 5 17 15 +17 1 tri3 13 17 7 +18 1 tri3 2 3 16 +19 1 tri3 14 4 15 +20 1 tri3 9 16 3 +21 1 tri3 15 17 14 + +# marked edges +# marker p1 p2 +1 6 8 +1 1 6 +1 2 11 +1 2 3 +1 4 3 +4 5 15 +1 0 7 +1 8 0 +1 11 1 +1 7 17 +4 15 4 +1 17 5 +", + ) + .unwrap(); + if SAVE_FIGURE { + let mut draw = Draw::new(); + draw.show_cells(true) + .show_point_ids(true) + .show_point_marker(true) + .show_cell_ids(true) + .show_cell_marker(true) + .show_edge_markers(true) + .set_view_flag(false) + .set_size(800.0, 800.0); + draw.all(&mesh, "/tmp/gemlab/test_search_features_2d_works_2.svg") + .unwrap(); + } + + let feat_bry = Features::new(&mesh, false); // no internal edges + let feat_int = Features::new(&mesh, true); // all edges including internal edges + + // cells_by_points ----------------------------------------------------------------------- + let cells_by_points = true; + + // points along x=1.0 (but the internal is not included) + let (point_ids, cell_ids, edges) = feat_bry + .search_features_2d(cells_by_points, At::X(1.0), |_| true) + .unwrap(); + assert_eq!(point_ids, vec![3, 4, 6]); + assert_eq!(cell_ids, vec![1, 2, 14, 18, 19]); // only boundary cells + assert_eq!(edges.all.iter().map(|e| e.key()).collect::>(), vec![(3, 4)]); + + // points along x=1.0 with y<=1.0 (but the internal is not included) + let (point_ids, cell_ids, edges) = feat_bry + .search_features_2d(cells_by_points, At::X(1.0), |x| x[1] <= 1.0) + .unwrap(); + assert_eq!(point_ids, vec![3, 6]); + assert_eq!(cell_ids, vec![1, 2, 14, 18]); + assert_eq!(edges.all.len(), 0); // no boundary edges + + // points along x=1.0 (with internal edges) + let (point_ids, cell_ids, edges) = feat_int + .search_features_2d(cells_by_points, At::X(1.0), |_| true) + .unwrap(); + assert_eq!(point_ids, vec![3, 4, 6, 9]); + assert_eq!(cell_ids, vec![1, 2, 4, 6, 7, 8, 9, 11, 14, 18, 19, 20]); + assert_eq!( + edges.all.iter().map(|e| e.key()).collect::>(), + vec![(3, 4), (3, 9), (6, 9)] + ); + + // points along x=1.0 with y<=1.0 (with internal edges) + let (point_ids, cell_ids, edges) = feat_int + .search_features_2d(cells_by_points, At::X(1.0), |x| x[1] <= 1.0) + .unwrap(); + assert_eq!(point_ids, vec![3, 6, 9]); + assert_eq!(cell_ids, vec![1, 2, 4, 6, 7, 8, 9, 11, 14, 18, 20]); + assert_eq!( + edges.all.iter().map(|e| e.key()).collect::>(), + vec![(3, 9), (6, 9)] + ); + + // cells_by_edges ----------------------------------------------------------------------- + let cells_by_points = false; + + // points along x=1.0 (but the internal is not included) + let (point_ids, cell_ids, edges) = feat_bry + .search_features_2d(cells_by_points, At::X(1.0), |_| true) + .unwrap(); + assert_eq!(point_ids, vec![3, 4, 6]); + assert_eq!(cell_ids, vec![14]); // only boundary cells attached to boundary edges + assert_eq!(edges.all.iter().map(|e| e.key()).collect::>(), vec![(3, 4)]); + + // points along x=1.0 with y<=1.0 (but the internal is not included) + let (point_ids, cell_ids, edges) = feat_bry + .search_features_2d(cells_by_points, At::X(1.0), |x| x[1] <= 1.0) + .unwrap(); + assert_eq!(point_ids, vec![3, 6]); + assert_eq!(cell_ids.len(), 0); // no boundary cells attached to boundary edges are found + assert_eq!(edges.all.len(), 0); // no boundary edges + + // points along x=1.0 (with internal edges) + let (point_ids, cell_ids, edges) = feat_int + .search_features_2d(cells_by_points, At::X(1.0), |_| true) + .unwrap(); + assert_eq!(point_ids, vec![3, 4, 6, 9]); + assert_eq!(cell_ids, vec![2, 6, 11, 14, 20]); + assert_eq!( + edges.all.iter().map(|e| e.key()).collect::>(), + vec![(3, 4), (3, 9), (6, 9)] + ); + + // points along x=1.0 with y<=1.0 (with internal edges) + let (point_ids, cell_ids, edges) = feat_int + .search_features_2d(cells_by_points, At::X(1.0), |x| x[1] <= 1.0) + .unwrap(); + assert_eq!(point_ids, vec![3, 6, 9]); + assert_eq!(cell_ids, vec![2, 6, 11, 20]); + assert_eq!( + edges.all.iter().map(|e| e.key()).collect::>(), + vec![(3, 9), (6, 9)] + ); } } diff --git a/src/mesh/grid_cells.rs b/src/mesh/grid_cells.rs index f67b487..cf6f638 100644 --- a/src/mesh/grid_cells.rs +++ b/src/mesh/grid_cells.rs @@ -638,10 +638,10 @@ mod tests { \n\ Histogram of container num items\n\ ================================\n\ - [0,1) | 0 \n\ - [1,2) | 0 \n\ - [2,3) | 4 ************************************************************\n\ - \x20\x20sum = 4\n" + [0, 1) | 0 \n\ + [1, 2) | 0 \n\ + [2, 3) | 4 ************************************************************\n\ + \x20\x20\x20sum = 4\n" ); let (mi, ma) = grid.limits(); assert_eq!(mi, xmin); diff --git a/src/mesh/mesh.rs b/src/mesh/mesh.rs index 9a75f15..d10f1d7 100644 --- a/src/mesh/mesh.rs +++ b/src/mesh/mesh.rs @@ -10,9 +10,6 @@ use std::fs::{self, File}; use std::io::BufReader; use std::path::Path; -/// Defines a tolerance to compare points in [Mesh::get_sorted_points()] -pub const TOL_COMPARE_POINTS: f64 = 1e-6; - /// Aliases usize as Point ID pub type PointId = usize; @@ -516,13 +513,6 @@ impl Mesh { /// # Notes /// /// 1. The filter is applied before sorting the points. - /// 2. The tolerance to compare points is [TOL_COMPARE_POINTS] times the range of the coordinates; i.e.: - /// - /// ```text - /// tol_x = TOL_COMPARE_POINTS * (xmax - xmin) - /// tol_y = TOL_COMPARE_POINTS * (ymax - ymin) - /// tol_z = TOL_COMPARE_POINTS * (zmax - zmin) - /// ``` /// /// # Output /// @@ -579,17 +569,10 @@ impl Mesh { } } } - // calculate tolerances to compare point coordinates - let (min, max) = self.get_limits(); - let mut tol = vec![TOL_COMPARE_POINTS; self.ndim]; - for i in 0..self.ndim { - tol[i] *= max[i] - min[i]; - } - // sort nodes by x โ†’ y โ†’ z let sorted_indices = if d3 { - argsort3_f64(&zz, &yy, &xx, &tol) + argsort3_f64(&zz, &yy, &xx) } else { - argsort2_f64(&yy, &xx, &tol) + argsort2_f64(&yy, &xx) }; sorted_indices.iter().map(|&i| k2id[&i]).collect() } diff --git a/src/recovery/get_extrap_matrix.rs b/src/recovery/get_extrap_matrix.rs index 7e9da1d..8b45356 100644 --- a/src/recovery/get_extrap_matrix.rs +++ b/src/recovery/get_extrap_matrix.rs @@ -823,7 +823,7 @@ mod tests { let u_nodal_original = Vector::from(&[1.0, 2.0, 3.0, 4.0]); // original U values at nodes let u_point = do_interpolate(&mut pad, &u_nodal_original, &gauss); // interpolated U values @ integration points let u_nodal = do_extrapolate(&ee, &u_point); - vec_approx_eq(&u_nodal, &u_nodal_original, 1e-15); + vec_approx_eq(&u_nodal, &u_nodal_original, 1e-14); } #[test] @@ -834,7 +834,7 @@ mod tests { let gauss = Gauss::new_sized(pad.kind.class(), 1).unwrap(); let ee = get_extrap_matrix(&mut pad, &gauss).unwrap(); let ee_correct = Matrix::from(&[[1.0], [1.0], [1.0], [1.0]]); // the pseudo-inverse of P = [0.25, 0.25, 0.25, 0.25] - mat_approx_eq(&ee, &ee_correct, 1e-15); + mat_approx_eq(&ee, &ee_correct, 1e-14); // check let u_nodal_original = Vector::from(&[1.0, 2.0, 3.0, 4.0]); // original U values at nodes diff --git a/zdocker/Dockerfile.Arch b/zdocker/Dockerfile.Arch new file mode 100644 index 0000000..250c865 --- /dev/null +++ b/zdocker/Dockerfile.Arch @@ -0,0 +1,21 @@ +FROM archlinux:base-devel + +# Update keyring first to avoid GPG errors on a fresh Arch image, +# then upgrade the system and install Rust and core dependencies. +RUN pacman -Sy --noconfirm archlinux-keyring && \ + pacman -Syu --noconfirm rust blas-openblas rust suitesparse + +# Create a non-root user for building and testing +RUN useradd -m user + +# Copy project files; --chown ensures the destination directory and all +# files are owned by user, so cargo can write target/ during the build +COPY --chown=user:user . /usr/local/src/gemlab/ +WORKDIR /usr/local/src/gemlab + +# Set the user +USER user + +# Test the code +ENV RUST_BACKTRACE=1 +RUN cargo test -- --nocapture diff --git a/zdocker/Dockerfile.Arch.Mkl.Local b/zdocker/Dockerfile.Arch.Mkl.Local new file mode 100644 index 0000000..c1cc979 --- /dev/null +++ b/zdocker/Dockerfile.Arch.Mkl.Local @@ -0,0 +1,48 @@ +FROM archlinux:base-devel + +# Update keyring first to avoid GPG errors on a fresh Arch image, +# then upgrade the system and install Rust and core dependencies; +# suitesparse is installed without blas-openblas because MKL will be used instead +RUN pacman -Sy --noconfirm archlinux-keyring && \ + pacman -Syu --noconfirm rust + +# Create a non-root user for building and testing +RUN useradd -m user + +# Give user password-less sudo so yay/makepkg can install AUR packages +RUN echo 'user ALL=(ALL) NOPASSWD: ALL' >> /etc/sudoers + +# Copy project files; --chown ensures the destination directory and all +# files are owned by user, so cargo can write target/ during the build +COPY --chown=user:user . /usr/local/src/gemlab/ +WORKDIR /usr/local/src/gemlab + +# Install the Intel oneAPI toolkit (MKL + compilers); +# must run as root because pacman needs it +RUN bash zscripts/arch-install-intel-toolkit.bash + +# Install metis from Chaotic-AUR (pre-built AUR packages); +# must run as root because pacman needs it +RUN bash -c "pacman-key --init && \ + pacman-key --recv-key 3056513887B78AEB --keyserver keyserver.ubuntu.com && \ + pacman-key --lsign-key 3056513887B78AEB && \ + pacman -U --noconfirm 'https://cdn-mirror.chaotic.cx/chaotic-aur/chaotic-keyring.pkg.tar.zst' && \ + pacman -U --noconfirm 'https://cdn-mirror.chaotic.cx/chaotic-aur/chaotic-mirrorlist.pkg.tar.zst' && \ + echo >> /etc/pacman.conf && \ + echo '[chaotic-aur]' >> /etc/pacman.conf && \ + echo 'Include = /etc/pacman.d/chaotic-mirrorlist' >> /etc/pacman.conf && \ + pacman -Sy --noconfirm metis" + +# Compile and install MUMPS with Intel MKL backend; +# must run as root because pacman and ldconfig need it +RUN bash zscripts/arch-compile-mumps.bash 1 + +# Compile and install SuiteSparse as well (with Intel MKL) +RUN bash zscripts/arch-compile-suitesparse.bash 1 + +# Set the user +USER user + +# MKL requires sourcing setvars.sh in the same shell as cargo test +ENV RUST_BACKTRACE=1 +RUN bash -c "source /opt/intel/oneapi/setvars.sh && cargo test --features intel_mkl,local_sparse" diff --git a/zscripts/arch-compile-mumps.bash b/zscripts/arch-compile-mumps.bash new file mode 100755 index 0000000..cf22220 --- /dev/null +++ b/zscripts/arch-compile-mumps.bash @@ -0,0 +1,102 @@ +#!/bin/bash + +# Exit on error (-e), treat unset variables as errors (-u), and propagate +# pipeline failures (-o pipefail) so any silent failure is caught early +set -euo pipefail + +# When running as root (e.g. inside Docker), act as a no-op wrapper; +# otherwise delegate to the real sudo +sudo () { + [[ $EUID = 0 ]] || set -- command sudo "$@" + "$@" +} + +# Optional first argument: set to 1 to build against Intel MKL/oneAPI, +# or 0 (default) to use the OpenBLAS + gfortran toolchain. +# When using MKL, run arch-install-intel-toolkit.bash first. +USE_INTEL_MKL=${1:-0} + +# Installation paths and source version +VERSION="5.9.0" +PREFIX="/usr/local" +INCDIR="$PREFIX/include/mumps" +LIBDIR="$PREFIX/lib/mumps" +PDIR="$(pwd)/zscripts/makefiles-mumps" + +# Install build tools and BLAS/LAPACK dependencies; +sudo pacman -S --noconfirm base-devel curl +if [ "${USE_INTEL_MKL}" = "1" ]; then + set +u # setvars.sh references variables before setting them; suppress -u temporarily + source /opt/intel/oneapi/setvars.sh # sets CC, FC, and PATH for the ifx compiler + set -u + export | grep -i MKLROOT +else + sudo pacman -S --noconfirm gcc-fortran openmp +fi + +# Metis is needed for graph-partitioning-based reordering; +# skip if already installed (e.g. via Chaotic-AUR in Docker) +if pacman -Q metis &>/dev/null; then + echo "... metis (package) already installed" +else + # yay must run as non-root; when this script is invoked as root (e.g. + # outside Docker), delegate to the unprivileged 'user' account that has + # NOPASSWD sudo + if [ $EUID = 0 ]; then + su - user -c "yay -S --noconfirm metis" + else + yay -S --noconfirm metis + fi +fi + +# Download the source tarball from Debian (reuse it if already present) +MUMPS_GZ="mumps_${VERSION}.orig.tar.gz" +MUMPS_DIR="MUMPS_${VERSION}" +cd /tmp +if [ -d "$MUMPS_DIR" ]; then + echo "... removing previous $MUMPS_DIR directory" + rm -rf "$MUMPS_DIR" +fi +if [ -f "$MUMPS_GZ" ]; then + echo "... using existing $MUMPS_GZ file" +else + curl -fL "http://deb.debian.org/debian/pool/main/m/mumps/${MUMPS_GZ}" -o "$MUMPS_GZ" +fi + +tar xzf "$MUMPS_GZ" +cd "$MUMPS_DIR" + +# Select the Makefile.inc that matches the chosen BLAS/LAPACK backend +if [ "${USE_INTEL_MKL}" = "1" ]; then + cp "$PDIR/MakefileMKL.inc" Makefile.inc +else + cp "$PDIR/Makefile.inc" Makefile.inc +fi + +# Create installation directories +sudo mkdir -p "$LIBDIR" "$INCDIR" + +# Build all four variants and copy the resulting libraries; +# make clean between targets is required because object files are not +# recompiled when only the output target changes +make d # double-precision, static +sudo cp -av lib/lib*.a "$LIBDIR/" + +make clean +make dshared # double-precision, shared +sudo cp -av lib/lib*.so "$LIBDIR/" + +make clean +make z # complex (double), static +sudo cp -av lib/lib*.a "$LIBDIR/" + +make clean +make zshared # complex (double), shared +sudo cp -av lib/lib*.so "$LIBDIR/" + +# Install public headers +sudo cp -av include/*.h "$INCDIR/" + +# Register the new library path with the dynamic linker +echo "$LIBDIR" | sudo tee /etc/ld.so.conf.d/mumps.conf >/dev/null +sudo ldconfig diff --git a/zscripts/arch-compile-suitesparse.bash b/zscripts/arch-compile-suitesparse.bash new file mode 100755 index 0000000..cba8685 --- /dev/null +++ b/zscripts/arch-compile-suitesparse.bash @@ -0,0 +1,82 @@ +#!/bin/bash + +# Exit on error (-e), treat unset variables as errors (-u), and propagate +# pipeline failures (-o pipefail) so any silent failure is caught early +set -euo pipefail + +# When running as root (e.g. inside Docker), act as a no-op wrapper; +# otherwise delegate to the real sudo +sudo () { + [[ $EUID = 0 ]] || set -- command sudo "$@" + "$@" +} + +# Optional first argument: set to 1 to build against Intel MKL/oneAPI, +# or 0 (default) to use the OpenBLAS toolchain. +# When using MKL, run arch-install-intel-toolkit.bash first. +USE_INTEL_MKL=${1:-0} + +# Installation paths +PREFIX="/usr/local" +INCDIR=$PREFIX/include/suitesparse +LIBDIR=$PREFIX/lib/suitesparse + +# Install build tools and BLAS/LAPACK dependencies +sudo pacman -S --noconfirm cmake curl git make +if [ "${USE_INTEL_MKL}" = "1" ]; then + set +u # setvars.sh references variables before setting them; suppress -u temporarily + source /opt/intel/oneapi/setvars.sh # sets CC, FC, and PATH for the ifx compiler + set -u + export | grep -i MKLROOT +else + sudo pacman -S --noconfirm clang openblas +fi + +# set the compiler/linker and cmake options for SuiteSparse +if [ "${USE_INTEL_MKL}" = "1" ]; then + COMP_LINK="-DCMAKE_C_COMPILER=icx -DCMAKE_CXX_COMPILER=icpx" + CMAKE_OPTIONS="${COMP_LINK} -DBLA_VENDOR=Intel10_64lp -DBLA_SIZEOF_INTEGER=4 -DSUITESPARSE_USE_FORTRAN=OFF" +else + COMP_LINK="-DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++" + CMAKE_OPTIONS="${COMP_LINK} -DBLA_VENDOR=OpenBLAS -DBLA_SIZEOF_INTEGER=4 -DSUITESPARSE_USE_FORTRAN=OFF" +fi + +# download the source code +cd /tmp +if ! [ -d "SuiteSparse" ]; then + git clone https://github.com/DrTimothyAldenDavis/SuiteSparse.git +fi +cd SuiteSparse + +# function to compile and copy results to local directories +action () { + local dir=$1 + cd $dir + make clean + CMAKE_OPTIONS=${CMAKE_OPTIONS} make local + CMAKE_OPTIONS=${CMAKE_OPTIONS} make install + cd .. +} + +# compile and copy results to local directories +action SuiteSparse_config +action AMD +action CAMD +action CCOLAMD +action COLAMD +action CHOLMOD +action UMFPACK +action BTF +action KLU + +# copy include files +sudo mkdir -p $INCDIR/ +sudo cp -av include/suitesparse/*.h $INCDIR/ + +# copy libray files +sudo mkdir -p $LIBDIR/ +sudo cp -av lib/* $LIBDIR/ + +# update ldconfig +echo "${LIBDIR}" | sudo tee /etc/ld.so.conf.d/suitesparse.conf >/dev/null +sudo ldconfig diff --git a/zscripts/arch-install-intel-toolkit.bash b/zscripts/arch-install-intel-toolkit.bash new file mode 100755 index 0000000..a1a833f --- /dev/null +++ b/zscripts/arch-install-intel-toolkit.bash @@ -0,0 +1,32 @@ +#!/bin/bash + +# Exit on error (-e), treat unset variables as errors (-u), and propagate +# pipeline failures (-o pipefail) so any silent failure is caught early +set -euo pipefail + +# When running as root (e.g. inside Docker), act as a no-op wrapper; +# otherwise delegate to the real sudo +sudo () { + [[ $EUID = 0 ]] || set -- command sudo "$@" + "$@" +} + +# Install the Intel oneAPI toolkit package +sudo pacman -S --noconfirm intel-oneapi-toolkit + +# Source the oneAPI environment to make MKLROOT and compiler paths available; +# suppress -u temporarily because setvars.sh references variables before setting them +set +u +source /opt/intel/oneapi/setvars.sh # sets CC, FC, MKLROOT, and PATH for the ifx compiler +set -u + +# Remove the broken non-ELF file that ships with MKL and confuses ldconfig +# See https://community.intel.com/t5/Intel-oneAPI-Math-Kernel-Library/ldconfig-opt-intel-oneapi-redist-lib-libmkl-sycl-so-is-not-an/m-p/1549240#M35528 +WEIRD_FILE="$MKLROOT/lib/intel64/libmkl_sycl.so" +if [ -f "$WEIRD_FILE" ]; then + sudo mv "$WEIRD_FILE" "$MKLROOT/lib/intel64/libmkl_sick.txt" +fi + +# Register the MKL library path with the dynamic linker +echo "$MKLROOT/lib/intel64" | sudo tee /etc/ld.so.conf.d/intel-mkl.conf >/dev/null +sudo ldconfig diff --git a/zscripts/docker-build-arch.bash b/zscripts/docker-build-arch.bash new file mode 100755 index 0000000..463342e --- /dev/null +++ b/zscripts/docker-build-arch.bash @@ -0,0 +1,24 @@ +#!/bin/bash + +set -euo pipefail + +ALTERNATIVE="${1:-0}" + +# image name and Dockerfile +if [ "${ALTERNATIVE}" = "1" ]; then + DOCKERFILE="zdocker/Dockerfile.Arch.Mkl.Local" + NAME="cpmech/gemlab_arch_mkl_local" +else + DOCKERFILE="zdocker/Dockerfile.Arch" + NAME="cpmech/gemlab_arch" +fi + +# build Docker image +docker build \ + -f "${DOCKERFILE}" \ + -t "${NAME}" \ + . + +echo +echo "... SUCCESS: image ${NAME} created ..." +echo diff --git a/zscripts/docker-shell-into.bash b/zscripts/docker-shell-into.bash new file mode 100755 index 0000000..e60b650 --- /dev/null +++ b/zscripts/docker-shell-into.bash @@ -0,0 +1,7 @@ +#!/bin/bash + +set -euo pipefail + +NAME="cpmech/gemlab_arch" + +docker run --rm -it "${NAME}:latest" /bin/bash diff --git a/zscripts/makefiles-mumps/MakefileMKL.inc b/zscripts/makefiles-mumps/MakefileMKL.inc new file mode 100644 index 0000000..5db78ec --- /dev/null +++ b/zscripts/makefiles-mumps/MakefileMKL.inc @@ -0,0 +1,70 @@ +# +# This file is part of MUMPS 5.7.1, released +# on Thu May 2 10:15:09 UTC 2024 +# + +# must be at the top +PLAT = _cpmech + +# Begin orderings +SCOTCHDIR = +ISCOTCH = + +LSCOTCH = + +LPORDDIR = $(topdir)/PORD/lib/ +IPORD = -I$(topdir)/PORD/include/ +LPORD = -L$(LPORDDIR) -lpord$(PLAT) + +LMETISDIR = /usr/lib +IMETIS = -I/usr/include + +LMETIS = -L$(LMETISDIR) -lmetis + +# Corresponding variables reused later +ORDERINGSF = -Dmetis -Dpord +ORDERINGSC = $(ORDERINGSF) + +LORDERINGS = $(LMETIS) $(LPORD) $(LSCOTCH) +IORDERINGSF = $(ISCOTCH) +IORDERINGSC = $(IMETIS) $(IPORD) $(ISCOTCH) +# End orderings +################################################################################ + +LIBEXT_SHARED = .so +SONAME = -soname +SHARED_OPT = -shared +FPIC_OPT = -fPIC +# Adapt/uncomment RPATH_OPT to avoid modifying +# LD_LIBRARY_PATH in case of shared libraries +# RPATH_OPT = -Wl,-rpath,/path/to/MUMPS_x.y.z/lib/ +LIBEXT = .a +OUTC = -o +OUTF = -o +RM = /bin/rm -f +CC = icx +FC = ifx +FL = ifx +AR = ar vr +RANLIB = ranlib +MKLROOT=/opt/intel/oneapi/mkl/latest/lib/intel64 +LAPACK = -L$(MKLROOT) -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core + +INCSEQ = -I$(topdir)/libseq +LIBSEQ = $(LAPACK) -L$(topdir)/libseq -lmpiseq$(PLAT) + +LIBBLAS = -L$(MKLROOT) -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core +LIBOTHERS = -lpthread + +#Preprocessor defs for calling Fortran from C (-DAdd_ or -DAdd__ or -DUPPER) +CDEFS = -DAdd_ + +#Begin Optimized options +OPTF = -O -nofor-main -qopenmp -DGEMMT_AVAILABLE +OPTL = -O -nofor-main -qopenmp +OPTC = -O -qopenmp +#End Optimized options + +INCS = $(INCSEQ) +LIBS = $(LIBSEQ) +LIBSEQNEEDED = libseqneeded