Thanks for using Compiler Explorer
Sponsors
Jakt
C++
Ada
Algol68
Analysis
Android Java
Android Kotlin
Assembly
C
C3
Carbon
C with Coccinelle
C++ with Coccinelle
C++ (Circle)
CIRCT
Clean
CMake
CMakeScript
COBOL
C++ for OpenCL
MLIR
Cppx
Cppx-Blue
Cppx-Gold
Cpp2-cppfront
Crystal
C#
CUDA C++
D
Dart
Elixir
Erlang
Fortran
F#
GLSL
Go
Haskell
HLSL
Hook
Hylo
IL
ispc
Java
Julia
Kotlin
LLVM IR
LLVM MIR
Modula-2
Mojo
Nim
Numba
Nix
Objective-C
Objective-C++
OCaml
Odin
OpenCL C
Pascal
Pony
PTX
Python
Racket
Raku
Ruby
Rust
Sail
Snowball
Scala
Slang
Solidity
Spice
SPIR-V
Swift
LLVM TableGen
Toit
Triton
TypeScript Native
V
Vala
Visual Basic
Vyper
WASM
Zig
Javascript
GIMPLE
Ygen
sway
c++ source #1
Output
Compile to binary object
Link to binary
Execute the code
Intel asm syntax
Demangle identifiers
Verbose demangling
Filters
Unused labels
Library functions
Directives
Comments
Horizontal whitespace
Debug intrinsics
Compiler
6502-c++ 11.1.0
ARM GCC 10.2.0
ARM GCC 10.3.0
ARM GCC 10.4.0
ARM GCC 10.5.0
ARM GCC 11.1.0
ARM GCC 11.2.0
ARM GCC 11.3.0
ARM GCC 11.4.0
ARM GCC 12.1.0
ARM GCC 12.2.0
ARM GCC 12.3.0
ARM GCC 12.4.0
ARM GCC 12.5.0
ARM GCC 13.1.0
ARM GCC 13.2.0
ARM GCC 13.2.0 (unknown-eabi)
ARM GCC 13.3.0
ARM GCC 13.3.0 (unknown-eabi)
ARM GCC 13.4.0
ARM GCC 13.4.0 (unknown-eabi)
ARM GCC 14.1.0
ARM GCC 14.1.0 (unknown-eabi)
ARM GCC 14.2.0
ARM GCC 14.2.0 (unknown-eabi)
ARM GCC 14.3.0
ARM GCC 14.3.0 (unknown-eabi)
ARM GCC 15.1.0
ARM GCC 15.1.0 (unknown-eabi)
ARM GCC 15.2.0
ARM GCC 15.2.0 (unknown-eabi)
ARM GCC 4.5.4
ARM GCC 4.6.4
ARM GCC 5.4
ARM GCC 6.3.0
ARM GCC 6.4.0
ARM GCC 7.3.0
ARM GCC 7.5.0
ARM GCC 8.2.0
ARM GCC 8.5.0
ARM GCC 9.3.0
ARM GCC 9.4.0
ARM GCC 9.5.0
ARM GCC trunk
ARM gcc 10.2.1 (none)
ARM gcc 10.3.1 (2021.07 none)
ARM gcc 10.3.1 (2021.10 none)
ARM gcc 11.2.1 (none)
ARM gcc 5.4.1 (none)
ARM gcc 7.2.1 (none)
ARM gcc 8.2 (WinCE)
ARM gcc 8.3.1 (none)
ARM gcc 9.2.1 (none)
ARM msvc v19.0 (ex-WINE)
ARM msvc v19.10 (ex-WINE)
ARM msvc v19.14 (ex-WINE)
ARM64 Morello gcc 10.1 Alpha 2
ARM64 gcc 10.2
ARM64 gcc 10.3
ARM64 gcc 10.4
ARM64 gcc 10.5.0
ARM64 gcc 11.1
ARM64 gcc 11.2
ARM64 gcc 11.3
ARM64 gcc 11.4.0
ARM64 gcc 12.1
ARM64 gcc 12.2.0
ARM64 gcc 12.3.0
ARM64 gcc 12.4.0
ARM64 gcc 12.5.0
ARM64 gcc 13.1.0
ARM64 gcc 13.2.0
ARM64 gcc 13.3.0
ARM64 gcc 13.4.0
ARM64 gcc 14.1.0
ARM64 gcc 14.2.0
ARM64 gcc 14.3.0
ARM64 gcc 15.1.0
ARM64 gcc 15.2.0
ARM64 gcc 4.9.4
ARM64 gcc 5.4
ARM64 gcc 5.5.0
ARM64 gcc 6.3
ARM64 gcc 6.4
ARM64 gcc 7.3
ARM64 gcc 7.5
ARM64 gcc 8.2
ARM64 gcc 8.5
ARM64 gcc 9.3
ARM64 gcc 9.4
ARM64 gcc 9.5
ARM64 gcc trunk
ARM64 msvc v19.14 (ex-WINE)
AVR gcc 10.3.0
AVR gcc 11.1.0
AVR gcc 12.1.0
AVR gcc 12.2.0
AVR gcc 12.3.0
AVR gcc 12.4.0
AVR gcc 12.5.0
AVR gcc 13.1.0
AVR gcc 13.2.0
AVR gcc 13.3.0
AVR gcc 13.4.0
AVR gcc 14.1.0
AVR gcc 14.2.0
AVR gcc 14.3.0
AVR gcc 15.1.0
AVR gcc 15.2.0
AVR gcc 4.5.4
AVR gcc 4.6.4
AVR gcc 5.4.0
AVR gcc 9.2.0
AVR gcc 9.3.0
Arduino Mega (1.8.9)
Arduino Uno (1.8.9)
BPF clang (trunk)
BPF clang 13.0.0
BPF clang 14.0.0
BPF clang 15.0.0
BPF clang 16.0.0
BPF clang 17.0.1
BPF clang 18.1.0
BPF clang 19.1.0
BPF clang 20.1.0
BPF clang 21.1.0
EDG (experimental reflection)
EDG 6.5
EDG 6.5 (GNU mode gcc 13)
EDG 6.6
EDG 6.6 (GNU mode gcc 13)
EDG 6.7
EDG 6.7 (GNU mode gcc 14)
FRC 2019
FRC 2020
FRC 2023
HPPA gcc 14.2.0
HPPA gcc 14.3.0
HPPA gcc 15.1.0
HPPA gcc 15.2.0
KVX ACB 4.1.0 (GCC 7.5.0)
KVX ACB 4.1.0-cd1 (GCC 7.5.0)
KVX ACB 4.10.0 (GCC 10.3.1)
KVX ACB 4.11.1 (GCC 10.3.1)
KVX ACB 4.12.0 (GCC 11.3.0)
KVX ACB 4.2.0 (GCC 7.5.0)
KVX ACB 4.3.0 (GCC 7.5.0)
KVX ACB 4.4.0 (GCC 7.5.0)
KVX ACB 4.6.0 (GCC 9.4.1)
KVX ACB 4.8.0 (GCC 9.4.1)
KVX ACB 4.9.0 (GCC 9.4.1)
KVX ACB 5.0.0 (GCC 12.2.1)
KVX ACB 5.2.0 (GCC 13.2.1)
LoongArch64 clang (trunk)
LoongArch64 clang 17.0.1
LoongArch64 clang 18.1.0
LoongArch64 clang 19.1.0
LoongArch64 clang 20.1.0
LoongArch64 clang 21.1.0
M68K gcc 13.1.0
M68K gcc 13.2.0
M68K gcc 13.3.0
M68K gcc 13.4.0
M68K gcc 14.1.0
M68K gcc 14.2.0
M68K gcc 14.3.0
M68K gcc 15.1.0
M68K gcc 15.2.0
M68k clang (trunk)
MRISC32 gcc (trunk)
MSP430 gcc 4.5.3
MSP430 gcc 5.3.0
MSP430 gcc 6.2.1
MinGW clang 14.0.3
MinGW clang 14.0.6
MinGW clang 15.0.7
MinGW clang 16.0.0
MinGW clang 16.0.2
MinGW gcc 11.3.0
MinGW gcc 12.1.0
MinGW gcc 12.2.0
MinGW gcc 13.1.0
MinGW gcc 14.3.0
MinGW gcc 15.2.0
RISC-V (32-bits) gcc (trunk)
RISC-V (32-bits) gcc 10.2.0
RISC-V (32-bits) gcc 10.3.0
RISC-V (32-bits) gcc 11.2.0
RISC-V (32-bits) gcc 11.3.0
RISC-V (32-bits) gcc 11.4.0
RISC-V (32-bits) gcc 12.1.0
RISC-V (32-bits) gcc 12.2.0
RISC-V (32-bits) gcc 12.3.0
RISC-V (32-bits) gcc 12.4.0
RISC-V (32-bits) gcc 12.5.0
RISC-V (32-bits) gcc 13.1.0
RISC-V (32-bits) gcc 13.2.0
RISC-V (32-bits) gcc 13.3.0
RISC-V (32-bits) gcc 13.4.0
RISC-V (32-bits) gcc 14.1.0
RISC-V (32-bits) gcc 14.2.0
RISC-V (32-bits) gcc 14.3.0
RISC-V (32-bits) gcc 15.1.0
RISC-V (32-bits) gcc 15.2.0
RISC-V (32-bits) gcc 8.2.0
RISC-V (32-bits) gcc 8.5.0
RISC-V (32-bits) gcc 9.4.0
RISC-V (64-bits) gcc (trunk)
RISC-V (64-bits) gcc 10.2.0
RISC-V (64-bits) gcc 10.3.0
RISC-V (64-bits) gcc 11.2.0
RISC-V (64-bits) gcc 11.3.0
RISC-V (64-bits) gcc 11.4.0
RISC-V (64-bits) gcc 12.1.0
RISC-V (64-bits) gcc 12.2.0
RISC-V (64-bits) gcc 12.3.0
RISC-V (64-bits) gcc 12.4.0
RISC-V (64-bits) gcc 12.5.0
RISC-V (64-bits) gcc 13.1.0
RISC-V (64-bits) gcc 13.2.0
RISC-V (64-bits) gcc 13.3.0
RISC-V (64-bits) gcc 13.4.0
RISC-V (64-bits) gcc 14.1.0
RISC-V (64-bits) gcc 14.2.0
RISC-V (64-bits) gcc 14.3.0
RISC-V (64-bits) gcc 15.1.0
RISC-V (64-bits) gcc 15.2.0
RISC-V (64-bits) gcc 8.2.0
RISC-V (64-bits) gcc 8.5.0
RISC-V (64-bits) gcc 9.4.0
RISC-V rv32gc clang (trunk)
RISC-V rv32gc clang 10.0.0
RISC-V rv32gc clang 10.0.1
RISC-V rv32gc clang 11.0.0
RISC-V rv32gc clang 11.0.1
RISC-V rv32gc clang 12.0.0
RISC-V rv32gc clang 12.0.1
RISC-V rv32gc clang 13.0.0
RISC-V rv32gc clang 13.0.1
RISC-V rv32gc clang 14.0.0
RISC-V rv32gc clang 15.0.0
RISC-V rv32gc clang 16.0.0
RISC-V rv32gc clang 17.0.1
RISC-V rv32gc clang 18.1.0
RISC-V rv32gc clang 19.1.0
RISC-V rv32gc clang 20.1.0
RISC-V rv32gc clang 21.1.0
RISC-V rv32gc clang 9.0.0
RISC-V rv32gc clang 9.0.1
RISC-V rv64gc clang (trunk)
RISC-V rv64gc clang 10.0.0
RISC-V rv64gc clang 10.0.1
RISC-V rv64gc clang 11.0.0
RISC-V rv64gc clang 11.0.1
RISC-V rv64gc clang 12.0.0
RISC-V rv64gc clang 12.0.1
RISC-V rv64gc clang 13.0.0
RISC-V rv64gc clang 13.0.1
RISC-V rv64gc clang 14.0.0
RISC-V rv64gc clang 15.0.0
RISC-V rv64gc clang 16.0.0
RISC-V rv64gc clang 17.0.1
RISC-V rv64gc clang 18.1.0
RISC-V rv64gc clang 19.1.0
RISC-V rv64gc clang 20.1.0
RISC-V rv64gc clang 21.1.0
RISC-V rv64gc clang 9.0.0
RISC-V rv64gc clang 9.0.1
Raspbian Buster
Raspbian Stretch
SPARC LEON gcc 12.2.0
SPARC LEON gcc 12.3.0
SPARC LEON gcc 12.4.0
SPARC LEON gcc 12.5.0
SPARC LEON gcc 13.1.0
SPARC LEON gcc 13.2.0
SPARC LEON gcc 13.3.0
SPARC LEON gcc 13.4.0
SPARC LEON gcc 14.1.0
SPARC LEON gcc 14.2.0
SPARC LEON gcc 14.3.0
SPARC LEON gcc 15.1.0
SPARC LEON gcc 15.2.0
SPARC gcc 12.2.0
SPARC gcc 12.3.0
SPARC gcc 12.4.0
SPARC gcc 12.5.0
SPARC gcc 13.1.0
SPARC gcc 13.2.0
SPARC gcc 13.3.0
SPARC gcc 13.4.0
SPARC gcc 14.1.0
SPARC gcc 14.2.0
SPARC gcc 14.3.0
SPARC gcc 15.1.0
SPARC gcc 15.2.0
SPARC64 gcc 12.2.0
SPARC64 gcc 12.3.0
SPARC64 gcc 12.4.0
SPARC64 gcc 12.5.0
SPARC64 gcc 13.1.0
SPARC64 gcc 13.2.0
SPARC64 gcc 13.3.0
SPARC64 gcc 13.4.0
SPARC64 gcc 14.1.0
SPARC64 gcc 14.2.0
SPARC64 gcc 14.3.0
SPARC64 gcc 15.1.0
SPARC64 gcc 15.2.0
TI C6x gcc 12.2.0
TI C6x gcc 12.3.0
TI C6x gcc 12.4.0
TI C6x gcc 12.5.0
TI C6x gcc 13.1.0
TI C6x gcc 13.2.0
TI C6x gcc 13.3.0
TI C6x gcc 13.4.0
TI C6x gcc 14.1.0
TI C6x gcc 14.2.0
TI C6x gcc 14.3.0
TI C6x gcc 15.1.0
TI C6x gcc 15.2.0
TI CL430 21.6.1
Tricore gcc 11.3.0 (EEESlab)
VAX gcc NetBSDELF 10.4.0
VAX gcc NetBSDELF 10.5.0 (Nov 15 03:50:22 2023)
VAX gcc NetBSDELF 12.4.0 (Apr 16 05:27 2025)
WebAssembly clang (trunk)
Xtensa ESP32 gcc 11.2.0 (2022r1)
Xtensa ESP32 gcc 12.2.0 (20230208)
Xtensa ESP32 gcc 14.2.0 (20241119)
Xtensa ESP32 gcc 8.2.0 (2019r2)
Xtensa ESP32 gcc 8.2.0 (2020r1)
Xtensa ESP32 gcc 8.2.0 (2020r2)
Xtensa ESP32 gcc 8.4.0 (2020r3)
Xtensa ESP32 gcc 8.4.0 (2021r1)
Xtensa ESP32 gcc 8.4.0 (2021r2)
Xtensa ESP32-S2 gcc 11.2.0 (2022r1)
Xtensa ESP32-S2 gcc 12.2.0 (20230208)
Xtensa ESP32-S2 gcc 14.2.0 (20241119)
Xtensa ESP32-S2 gcc 8.2.0 (2019r2)
Xtensa ESP32-S2 gcc 8.2.0 (2020r1)
Xtensa ESP32-S2 gcc 8.2.0 (2020r2)
Xtensa ESP32-S2 gcc 8.4.0 (2020r3)
Xtensa ESP32-S2 gcc 8.4.0 (2021r1)
Xtensa ESP32-S2 gcc 8.4.0 (2021r2)
Xtensa ESP32-S3 gcc 11.2.0 (2022r1)
Xtensa ESP32-S3 gcc 12.2.0 (20230208)
Xtensa ESP32-S3 gcc 14.2.0 (20241119)
Xtensa ESP32-S3 gcc 8.4.0 (2020r3)
Xtensa ESP32-S3 gcc 8.4.0 (2021r1)
Xtensa ESP32-S3 gcc 8.4.0 (2021r2)
arm64 msvc v19.20 VS16.0
arm64 msvc v19.21 VS16.1
arm64 msvc v19.22 VS16.2
arm64 msvc v19.23 VS16.3
arm64 msvc v19.24 VS16.4
arm64 msvc v19.25 VS16.5
arm64 msvc v19.27 VS16.7
arm64 msvc v19.28 VS16.8
arm64 msvc v19.28 VS16.9
arm64 msvc v19.29 VS16.10
arm64 msvc v19.29 VS16.11
arm64 msvc v19.30 VS17.0
arm64 msvc v19.31 VS17.1
arm64 msvc v19.32 VS17.2
arm64 msvc v19.33 VS17.3
arm64 msvc v19.34 VS17.4
arm64 msvc v19.35 VS17.5
arm64 msvc v19.36 VS17.6
arm64 msvc v19.37 VS17.7
arm64 msvc v19.38 VS17.8
arm64 msvc v19.39 VS17.9
arm64 msvc v19.40 VS17.10
arm64 msvc v19.41 VS17.11
arm64 msvc v19.42 VS17.12
arm64 msvc v19.43 VS17.13
arm64 msvc v19.latest
armv7-a clang (trunk)
armv7-a clang 10.0.0
armv7-a clang 10.0.1
armv7-a clang 11.0.0
armv7-a clang 11.0.1
armv7-a clang 12.0.0
armv7-a clang 12.0.1
armv7-a clang 13.0.0
armv7-a clang 13.0.1
armv7-a clang 14.0.0
armv7-a clang 15.0.0
armv7-a clang 16.0.0
armv7-a clang 17.0.1
armv7-a clang 18.1.0
armv7-a clang 19.1.0
armv7-a clang 20.1.0
armv7-a clang 21.1.0
armv7-a clang 9.0.0
armv7-a clang 9.0.1
armv8-a clang (all architectural features, trunk)
armv8-a clang (trunk)
armv8-a clang 10.0.0
armv8-a clang 10.0.1
armv8-a clang 11.0.0
armv8-a clang 11.0.1
armv8-a clang 12.0.0
armv8-a clang 13.0.0
armv8-a clang 14.0.0
armv8-a clang 15.0.0
armv8-a clang 16.0.0
armv8-a clang 17.0.1
armv8-a clang 18.1.0
armv8-a clang 19.1.0
armv8-a clang 20.1.0
armv8-a clang 21.1.0
armv8-a clang 9.0.0
armv8-a clang 9.0.1
clad trunk (clang 21.1.0)
clad v1.10 (clang 20.1.0)
clad v1.8 (clang 18.1.0)
clad v1.9 (clang 19.1.0)
clad v2.00 (clang 20.1.0)
clang-cl 18.1.0
ellcc 0.1.33
ellcc 0.1.34
ellcc 2017-07-16
ez80-clang 15.0.0
ez80-clang 15.0.7
hexagon-clang 16.0.5
llvm-mos atari2600-3e
llvm-mos atari2600-4k
llvm-mos atari2600-common
llvm-mos atari5200-supercart
llvm-mos atari8-cart-megacart
llvm-mos atari8-cart-std
llvm-mos atari8-cart-xegs
llvm-mos atari8-common
llvm-mos atari8-dos
llvm-mos c128
llvm-mos c64
llvm-mos commodore
llvm-mos cpm65
llvm-mos cx16
llvm-mos dodo
llvm-mos eater
llvm-mos mega65
llvm-mos nes
llvm-mos nes-action53
llvm-mos nes-cnrom
llvm-mos nes-gtrom
llvm-mos nes-mmc1
llvm-mos nes-mmc3
llvm-mos nes-nrom
llvm-mos nes-unrom
llvm-mos nes-unrom-512
llvm-mos osi-c1p
llvm-mos pce
llvm-mos pce-cd
llvm-mos pce-common
llvm-mos pet
llvm-mos rp6502
llvm-mos rpc8e
llvm-mos supervision
llvm-mos vic20
loongarch64 gcc 12.2.0
loongarch64 gcc 12.3.0
loongarch64 gcc 12.4.0
loongarch64 gcc 12.5.0
loongarch64 gcc 13.1.0
loongarch64 gcc 13.2.0
loongarch64 gcc 13.3.0
loongarch64 gcc 13.4.0
loongarch64 gcc 14.1.0
loongarch64 gcc 14.2.0
loongarch64 gcc 14.3.0
loongarch64 gcc 15.1.0
loongarch64 gcc 15.2.0
mips clang 13.0.0
mips clang 14.0.0
mips clang 15.0.0
mips clang 16.0.0
mips clang 17.0.1
mips clang 18.1.0
mips clang 19.1.0
mips clang 20.1.0
mips clang 21.1.0
mips gcc 11.2.0
mips gcc 12.1.0
mips gcc 12.2.0
mips gcc 12.3.0
mips gcc 12.4.0
mips gcc 12.5.0
mips gcc 13.1.0
mips gcc 13.2.0
mips gcc 13.3.0
mips gcc 13.4.0
mips gcc 14.1.0
mips gcc 14.2.0
mips gcc 14.3.0
mips gcc 15.1.0
mips gcc 15.2.0
mips gcc 4.9.4
mips gcc 5.4
mips gcc 5.5.0
mips gcc 9.3.0 (codescape)
mips gcc 9.5.0
mips64 (el) gcc 12.1.0
mips64 (el) gcc 12.2.0
mips64 (el) gcc 12.3.0
mips64 (el) gcc 12.4.0
mips64 (el) gcc 12.5.0
mips64 (el) gcc 13.1.0
mips64 (el) gcc 13.2.0
mips64 (el) gcc 13.3.0
mips64 (el) gcc 13.4.0
mips64 (el) gcc 14.1.0
mips64 (el) gcc 14.2.0
mips64 (el) gcc 14.3.0
mips64 (el) gcc 15.1.0
mips64 (el) gcc 15.2.0
mips64 (el) gcc 4.9.4
mips64 (el) gcc 5.4.0
mips64 (el) gcc 5.5.0
mips64 (el) gcc 9.5.0
mips64 clang 13.0.0
mips64 clang 14.0.0
mips64 clang 15.0.0
mips64 clang 16.0.0
mips64 clang 17.0.1
mips64 clang 18.1.0
mips64 clang 19.1.0
mips64 clang 20.1.0
mips64 clang 21.1.0
mips64 gcc 11.2.0
mips64 gcc 12.1.0
mips64 gcc 12.2.0
mips64 gcc 12.3.0
mips64 gcc 12.4.0
mips64 gcc 12.5.0
mips64 gcc 13.1.0
mips64 gcc 13.2.0
mips64 gcc 13.3.0
mips64 gcc 13.4.0
mips64 gcc 14.1.0
mips64 gcc 14.2.0
mips64 gcc 14.3.0
mips64 gcc 15.1.0
mips64 gcc 15.2.0
mips64 gcc 4.9.4
mips64 gcc 5.4.0
mips64 gcc 5.5.0
mips64 gcc 9.5.0
mips64el clang 13.0.0
mips64el clang 14.0.0
mips64el clang 15.0.0
mips64el clang 16.0.0
mips64el clang 17.0.1
mips64el clang 18.1.0
mips64el clang 19.1.0
mips64el clang 20.1.0
mips64el clang 21.1.0
mipsel clang 13.0.0
mipsel clang 14.0.0
mipsel clang 15.0.0
mipsel clang 16.0.0
mipsel clang 17.0.1
mipsel clang 18.1.0
mipsel clang 19.1.0
mipsel clang 20.1.0
mipsel clang 21.1.0
mipsel gcc 12.1.0
mipsel gcc 12.2.0
mipsel gcc 12.3.0
mipsel gcc 12.4.0
mipsel gcc 12.5.0
mipsel gcc 13.1.0
mipsel gcc 13.2.0
mipsel gcc 13.3.0
mipsel gcc 13.4.0
mipsel gcc 14.1.0
mipsel gcc 14.2.0
mipsel gcc 14.3.0
mipsel gcc 15.1.0
mipsel gcc 15.2.0
mipsel gcc 4.9.4
mipsel gcc 5.4.0
mipsel gcc 5.5.0
mipsel gcc 9.5.0
nanoMIPS gcc 6.3.0 (mtk)
power gcc 11.2.0
power gcc 12.1.0
power gcc 12.2.0
power gcc 12.3.0
power gcc 12.4.0
power gcc 12.5.0
power gcc 13.1.0
power gcc 13.2.0
power gcc 13.3.0
power gcc 13.4.0
power gcc 14.1.0
power gcc 14.2.0
power gcc 14.3.0
power gcc 15.1.0
power gcc 15.2.0
power gcc 4.8.5
power64 AT12.0 (gcc8)
power64 AT13.0 (gcc9)
power64 gcc 11.2.0
power64 gcc 12.1.0
power64 gcc 12.2.0
power64 gcc 12.3.0
power64 gcc 12.4.0
power64 gcc 12.5.0
power64 gcc 13.1.0
power64 gcc 13.2.0
power64 gcc 13.3.0
power64 gcc 13.4.0
power64 gcc 14.1.0
power64 gcc 14.2.0
power64 gcc 14.3.0
power64 gcc 15.1.0
power64 gcc 15.2.0
power64 gcc trunk
power64le AT12.0 (gcc8)
power64le AT13.0 (gcc9)
power64le clang (trunk)
power64le gcc 11.2.0
power64le gcc 12.1.0
power64le gcc 12.2.0
power64le gcc 12.3.0
power64le gcc 12.4.0
power64le gcc 12.5.0
power64le gcc 13.1.0
power64le gcc 13.2.0
power64le gcc 13.3.0
power64le gcc 13.4.0
power64le gcc 14.1.0
power64le gcc 14.2.0
power64le gcc 14.3.0
power64le gcc 15.1.0
power64le gcc 15.2.0
power64le gcc 6.3.0
power64le gcc trunk
powerpc64 clang (trunk)
qnx 8.0.0
s390x gcc 11.2.0
s390x gcc 12.1.0
s390x gcc 12.2.0
s390x gcc 12.3.0
s390x gcc 12.4.0
s390x gcc 12.5.0
s390x gcc 13.1.0
s390x gcc 13.2.0
s390x gcc 13.3.0
s390x gcc 13.4.0
s390x gcc 14.1.0
s390x gcc 14.2.0
s390x gcc 14.3.0
s390x gcc 15.1.0
s390x gcc 15.2.0
sh gcc 12.2.0
sh gcc 12.3.0
sh gcc 12.4.0
sh gcc 12.5.0
sh gcc 13.1.0
sh gcc 13.2.0
sh gcc 13.3.0
sh gcc 13.4.0
sh gcc 14.1.0
sh gcc 14.2.0
sh gcc 14.3.0
sh gcc 15.1.0
sh gcc 15.2.0
sh gcc 4.9.4
sh gcc 9.5.0
vast (trunk)
x64 msvc v19.0 (ex-WINE)
x64 msvc v19.10 (ex-WINE)
x64 msvc v19.14 (ex-WINE)
x64 msvc v19.20 VS16.0
x64 msvc v19.21 VS16.1
x64 msvc v19.22 VS16.2
x64 msvc v19.23 VS16.3
x64 msvc v19.24 VS16.4
x64 msvc v19.25 VS16.5
x64 msvc v19.27 VS16.7
x64 msvc v19.28 VS16.8
x64 msvc v19.28 VS16.9
x64 msvc v19.29 VS16.10
x64 msvc v19.29 VS16.11
x64 msvc v19.30 VS17.0
x64 msvc v19.31 VS17.1
x64 msvc v19.32 VS17.2
x64 msvc v19.33 VS17.3
x64 msvc v19.34 VS17.4
x64 msvc v19.35 VS17.5
x64 msvc v19.36 VS17.6
x64 msvc v19.37 VS17.7
x64 msvc v19.38 VS17.8
x64 msvc v19.39 VS17.9
x64 msvc v19.40 VS17.10
x64 msvc v19.41 VS17.11
x64 msvc v19.42 VS17.12
x64 msvc v19.43 VS17.13
x64 msvc v19.latest
x86 djgpp 4.9.4
x86 djgpp 5.5.0
x86 djgpp 6.4.0
x86 djgpp 7.2.0
x86 msvc v19.0 (ex-WINE)
x86 msvc v19.10 (ex-WINE)
x86 msvc v19.14 (ex-WINE)
x86 msvc v19.20 VS16.0
x86 msvc v19.21 VS16.1
x86 msvc v19.22 VS16.2
x86 msvc v19.23 VS16.3
x86 msvc v19.24 VS16.4
x86 msvc v19.25 VS16.5
x86 msvc v19.27 VS16.7
x86 msvc v19.28 VS16.8
x86 msvc v19.28 VS16.9
x86 msvc v19.29 VS16.10
x86 msvc v19.29 VS16.11
x86 msvc v19.30 VS17.0
x86 msvc v19.31 VS17.1
x86 msvc v19.32 VS17.2
x86 msvc v19.33 VS17.3
x86 msvc v19.34 VS17.4
x86 msvc v19.35 VS17.5
x86 msvc v19.36 VS17.6
x86 msvc v19.37 VS17.7
x86 msvc v19.38 VS17.8
x86 msvc v19.39 VS17.9
x86 msvc v19.40 VS17.10
x86 msvc v19.41 VS17.11
x86 msvc v19.42 VS17.12
x86 msvc v19.43 VS17.13
x86 msvc v19.latest
x86 nvc++ 22.11
x86 nvc++ 22.7
x86 nvc++ 22.9
x86 nvc++ 23.1
x86 nvc++ 23.11
x86 nvc++ 23.3
x86 nvc++ 23.5
x86 nvc++ 23.7
x86 nvc++ 23.9
x86 nvc++ 24.1
x86 nvc++ 24.11
x86 nvc++ 24.3
x86 nvc++ 24.5
x86 nvc++ 24.7
x86 nvc++ 24.9
x86 nvc++ 25.1
x86 nvc++ 25.3
x86 nvc++ 25.5
x86 nvc++ 25.7
x86-64 Zapcc 190308
x86-64 clang (-fimplicit-constexpr)
x86-64 clang (Chris Bazley N3089)
x86-64 clang (EricWF contracts)
x86-64 clang (amd-staging)
x86-64 clang (assertions trunk)
x86-64 clang (clangir)
x86-64 clang (experimental -Wlifetime)
x86-64 clang (experimental P1061)
x86-64 clang (experimental P1144)
x86-64 clang (experimental P1221)
x86-64 clang (experimental P2998)
x86-64 clang (experimental P3068)
x86-64 clang (experimental P3309)
x86-64 clang (experimental P3367)
x86-64 clang (experimental P3372)
x86-64 clang (experimental P3385)
x86-64 clang (experimental P3776)
x86-64 clang (experimental metaprogramming - P2632)
x86-64 clang (old concepts branch)
x86-64 clang (p1974)
x86-64 clang (pattern matching - P2688)
x86-64 clang (reflection - C++26)
x86-64 clang (reflection - TS)
x86-64 clang (resugar)
x86-64 clang (string interpolation - P3412)
x86-64 clang (thephd.dev)
x86-64 clang (trunk)
x86-64 clang (variadic friends - P2893)
x86-64 clang (widberg)
x86-64 clang 10.0.0
x86-64 clang 10.0.0 (assertions)
x86-64 clang 10.0.1
x86-64 clang 11.0.0
x86-64 clang 11.0.0 (assertions)
x86-64 clang 11.0.1
x86-64 clang 12.0.0
x86-64 clang 12.0.0 (assertions)
x86-64 clang 12.0.1
x86-64 clang 13.0.0
x86-64 clang 13.0.0 (assertions)
x86-64 clang 13.0.1
x86-64 clang 14.0.0
x86-64 clang 14.0.0 (assertions)
x86-64 clang 15.0.0
x86-64 clang 15.0.0 (assertions)
x86-64 clang 16.0.0
x86-64 clang 16.0.0 (assertions)
x86-64 clang 17.0.1
x86-64 clang 17.0.1 (assertions)
x86-64 clang 18.1.0
x86-64 clang 18.1.0 (assertions)
x86-64 clang 19.1.0
x86-64 clang 19.1.0 (assertions)
x86-64 clang 2.6.0 (assertions)
x86-64 clang 2.7.0 (assertions)
x86-64 clang 2.8.0 (assertions)
x86-64 clang 2.9.0 (assertions)
x86-64 clang 20.1.0
x86-64 clang 20.1.0 (assertions)
x86-64 clang 21.1.0
x86-64 clang 21.1.0 (assertions)
x86-64 clang 3.0.0
x86-64 clang 3.0.0 (assertions)
x86-64 clang 3.1
x86-64 clang 3.1 (assertions)
x86-64 clang 3.2
x86-64 clang 3.2 (assertions)
x86-64 clang 3.3
x86-64 clang 3.3 (assertions)
x86-64 clang 3.4 (assertions)
x86-64 clang 3.4.1
x86-64 clang 3.5
x86-64 clang 3.5 (assertions)
x86-64 clang 3.5.1
x86-64 clang 3.5.2
x86-64 clang 3.6
x86-64 clang 3.6 (assertions)
x86-64 clang 3.7
x86-64 clang 3.7 (assertions)
x86-64 clang 3.7.1
x86-64 clang 3.8
x86-64 clang 3.8 (assertions)
x86-64 clang 3.8.1
x86-64 clang 3.9.0
x86-64 clang 3.9.0 (assertions)
x86-64 clang 3.9.1
x86-64 clang 4.0.0
x86-64 clang 4.0.0 (assertions)
x86-64 clang 4.0.1
x86-64 clang 5.0.0
x86-64 clang 5.0.0 (assertions)
x86-64 clang 5.0.1
x86-64 clang 5.0.2
x86-64 clang 6.0.0
x86-64 clang 6.0.0 (assertions)
x86-64 clang 6.0.1
x86-64 clang 7.0.0
x86-64 clang 7.0.0 (assertions)
x86-64 clang 7.0.1
x86-64 clang 7.1.0
x86-64 clang 8.0.0
x86-64 clang 8.0.0 (assertions)
x86-64 clang 8.0.1
x86-64 clang 9.0.0
x86-64 clang 9.0.0 (assertions)
x86-64 clang 9.0.1
x86-64 clang rocm-4.5.2
x86-64 clang rocm-5.0.2
x86-64 clang rocm-5.1.3
x86-64 clang rocm-5.2.3
x86-64 clang rocm-5.3.3
x86-64 clang rocm-5.7.0
x86-64 clang rocm-6.0.2
x86-64 clang rocm-6.1.2
x86-64 clang rocm-6.2.4
x86-64 clang rocm-6.3.3
x86-64 clang rocm-6.4.0
x86-64 gcc (P2034 lambdas)
x86-64 gcc (contract labels)
x86-64 gcc (contracts natural syntax)
x86-64 gcc (contracts)
x86-64 gcc (coroutines)
x86-64 gcc (modules)
x86-64 gcc (trunk)
x86-64 gcc 10.1
x86-64 gcc 10.2
x86-64 gcc 10.3
x86-64 gcc 10.3 (assertions)
x86-64 gcc 10.4
x86-64 gcc 10.4 (assertions)
x86-64 gcc 10.5
x86-64 gcc 10.5 (assertions)
x86-64 gcc 11.1
x86-64 gcc 11.1 (assertions)
x86-64 gcc 11.2
x86-64 gcc 11.2 (assertions)
x86-64 gcc 11.3
x86-64 gcc 11.3 (assertions)
x86-64 gcc 11.4
x86-64 gcc 11.4 (assertions)
x86-64 gcc 12.1
x86-64 gcc 12.1 (assertions)
x86-64 gcc 12.2
x86-64 gcc 12.2 (assertions)
x86-64 gcc 12.3
x86-64 gcc 12.3 (assertions)
x86-64 gcc 12.4
x86-64 gcc 12.4 (assertions)
x86-64 gcc 12.5
x86-64 gcc 12.5 (assertions)
x86-64 gcc 13.1
x86-64 gcc 13.1 (assertions)
x86-64 gcc 13.2
x86-64 gcc 13.2 (assertions)
x86-64 gcc 13.3
x86-64 gcc 13.3 (assertions)
x86-64 gcc 13.4
x86-64 gcc 13.4 (assertions)
x86-64 gcc 14.1
x86-64 gcc 14.1 (assertions)
x86-64 gcc 14.2
x86-64 gcc 14.2 (assertions)
x86-64 gcc 14.3
x86-64 gcc 14.3 (assertions)
x86-64 gcc 15.1
x86-64 gcc 15.1 (assertions)
x86-64 gcc 15.2
x86-64 gcc 15.2 (assertions)
x86-64 gcc 3.4.6
x86-64 gcc 4.0.4
x86-64 gcc 4.1.2
x86-64 gcc 4.4.7
x86-64 gcc 4.5.3
x86-64 gcc 4.6.4
x86-64 gcc 4.7.1
x86-64 gcc 4.7.2
x86-64 gcc 4.7.3
x86-64 gcc 4.7.4
x86-64 gcc 4.8.1
x86-64 gcc 4.8.2
x86-64 gcc 4.8.3
x86-64 gcc 4.8.4
x86-64 gcc 4.8.5
x86-64 gcc 4.9.0
x86-64 gcc 4.9.1
x86-64 gcc 4.9.2
x86-64 gcc 4.9.3
x86-64 gcc 4.9.4
x86-64 gcc 5.1
x86-64 gcc 5.2
x86-64 gcc 5.3
x86-64 gcc 5.4
x86-64 gcc 5.5
x86-64 gcc 6.1
x86-64 gcc 6.2
x86-64 gcc 6.3
x86-64 gcc 6.4
x86-64 gcc 6.5
x86-64 gcc 7.1
x86-64 gcc 7.2
x86-64 gcc 7.3
x86-64 gcc 7.4
x86-64 gcc 7.5
x86-64 gcc 8.1
x86-64 gcc 8.2
x86-64 gcc 8.3
x86-64 gcc 8.4
x86-64 gcc 8.5
x86-64 gcc 9.1
x86-64 gcc 9.2
x86-64 gcc 9.3
x86-64 gcc 9.4
x86-64 gcc 9.5
x86-64 icc 13.0.1
x86-64 icc 16.0.3
x86-64 icc 17.0.0
x86-64 icc 18.0.0
x86-64 icc 19.0.0
x86-64 icc 19.0.1
x86-64 icc 2021.1.2
x86-64 icc 2021.10.0
x86-64 icc 2021.2.0
x86-64 icc 2021.3.0
x86-64 icc 2021.4.0
x86-64 icc 2021.5.0
x86-64 icc 2021.6.0
x86-64 icc 2021.7.0
x86-64 icc 2021.7.1
x86-64 icc 2021.8.0
x86-64 icc 2021.9.0
x86-64 icx 2021.1.2
x86-64 icx 2021.2.0
x86-64 icx 2021.3.0
x86-64 icx 2021.4.0
x86-64 icx 2022.0.0
x86-64 icx 2022.1.0
x86-64 icx 2022.2.0
x86-64 icx 2022.2.1
x86-64 icx 2023.0.0
x86-64 icx 2023.1.0
x86-64 icx 2023.2.1
x86-64 icx 2024.0.0
x86-64 icx 2024.1.0
x86-64 icx 2024.2.0
x86-64 icx 2024.2.1
x86-64 icx 2025.0.0
x86-64 icx 2025.0.1
x86-64 icx 2025.0.3
x86-64 icx 2025.0.4
x86-64 icx 2025.1.0
x86-64 icx 2025.1.1
x86-64 icx 2025.2.0
x86-64 icx 2025.2.1
x86-64 icx 2025.2.1
z180-clang 15.0.0
z180-clang 15.0.7
z80-clang 15.0.0
z80-clang 15.0.7
zig c++ 0.10.0
zig c++ 0.11.0
zig c++ 0.12.0
zig c++ 0.12.1
zig c++ 0.13.0
zig c++ 0.14.0
zig c++ 0.14.1
zig c++ 0.15.1
zig c++ 0.6.0
zig c++ 0.7.0
zig c++ 0.7.1
zig c++ 0.8.0
zig c++ 0.9.0
zig c++ trunk
Options
Source code
#include <array> #include <bit> #include <cstddef> #include <cmath> #include <iostream> #include <utility> #include <valarray> #include <type_traits> #include <concepts> #include <string> #include <string_view> #include <vector> #include <map> #include <memory> #include <limits> #include <functional> #include <type_traits> #include <range/v3/view/zip.hpp> template<typename T, typename U> auto inline constexpr lerp(T const &lhs, U const &rhs, double const f) { return (1.0 - f) * lhs + f * rhs; } template<std::unsigned_integral auto n, typename T> requires (n > 0u) decltype(auto) inline constexpr pow(T const &x) { if constexpr (n == 1u) return x; else return x * pow<n - 1u>(x); } namespace VF { template<std::size_t d, std::size_t r1> struct TTensor { TTensor<d, r1 - 1u> TensorArr[d]; private: template<std::size_t... i> TTensor<d, r1> static constexpr Identidad(std::index_sequence<i...>) { return {(i / d == i % d ? 1.0 : 0.0)...}; } template<std::size_t... i> TTensor<d, r1> constexpr T(std::index_sequence<i...>) const { return {Slice<r1 - 1u, d, i>()...}; } template<std::size_t r2, std::size_t Stride, std::size_t Offset, std::size_t... i> TTensor<d, r2> constexpr Slice(std::index_sequence<i...>) const { return {std::bit_cast<std::array<double, pow<r1>(d)>>(*this)[i * Stride + Offset]...}; } template<std::size_t r2, std::size_t Stride, std::size_t Offset> TTensor<d, r2> constexpr Slice() const { return Slice<r2, Stride, Offset>(std::make_index_sequence<pow<r2>(d)>{}); } template<std::size_t... i> TTensor<d, r1> constexpr Suma(std::index_sequence<i...>, TTensor<d, r1> const &Tensor) const { return {(TensorArr[i] + Tensor[i])...}; } template<std::size_t... i> void constexpr SumaAsigna(std::index_sequence<i...>, TTensor<d, r1> const &Tensor) { ((TensorArr[i] += Tensor[i]), ...); } template<std::size_t... i> TTensor<d, r1> constexpr Resta(std::index_sequence<i...>, TTensor<d, r1> const &Tensor) const { return {(TensorArr[i] - Tensor[i])...}; } template<std::size_t... i> void constexpr RestaAsigna(std::index_sequence<i...>, TTensor<d, r1> const &Tensor) { ((TensorArr[i] -= Tensor[i]), ...); } template<std::size_t... i> TTensor<d, r1> constexpr Negativo(std::index_sequence<i...>) const { return {(-TensorArr[i])...}; } template<std::size_t r2, std::size_t... i> TTensor<d, r1 + r2> constexpr ProdExt(std::index_sequence<i...>, TTensor<d, r2> const &Tensor) const { return {(TensorArr[i] * Tensor)...}; } template<std::size_t... i> void constexpr ProdAsigna(std::index_sequence<i...>, double const Escalar) { ((TensorArr[i] *= Escalar), ...); } template<std::size_t... i> void constexpr DivAsigna(std::index_sequence<i...>, double const Escalar) { ((TensorArr[i] /= Escalar), ...); } template<std::size_t r2, std::size_t... i> TTensor<d, r1 + r2 - 2u> constexpr ProdInt(std::index_sequence<i...>, TTensor<d, r2> const &Tensor) const { return (... + (Slice<r1 - 1u, d, i>() * Tensor[i])); } template<std::size_t r2, std::size_t... i> requires (r1 == 1u) TTensor<d, r1 + r2 - 2u> constexpr ProdInt(std::index_sequence<i...>, TTensor<d, r2> const &Tensor) const { return (... + (TensorArr[i] * Tensor[i])); } template<std::size_t r2, std::size_t... i> TTensor<d, r1 + r2 - 4u> constexpr ProdIntDbl(std::index_sequence<i...>, TTensor<d, r2> const &Tensor) const { return (... + (Slice<r1 - 2u, d * d, i>() * Tensor.template Slice<r2 - 2u, 1u, i>())); } template<std::size_t... i> requires (r1 == 2u) TTensor<d, r1 - 2u> constexpr ProdIntDbl(std::index_sequence<i...>, TTensor<d, 2u> const &Tensor) const { return (... + (TensorArr[i] & Tensor[i])); } template<std::size_t... i> std::istream &ReadASCII(std::index_sequence<i...>, std::istream &is) { return (is >> ... >> TensorArr[i]); } template<std::size_t... i> requires (r1 == 1u) std::ostream &WriteASCII(std::index_sequence<i...>, std::ostream &os) const { return ((os << TensorArr[i] << ' '), ...) << TensorArr[d - 1u] << std::endl; } template<std::size_t... i> std::ostream &WriteASCII(std::index_sequence<i...>, std::ostream &os) const { return (os << ... << TensorArr[i]) << TensorArr[d - 1u] << std::endl; } public: TTensor<d, r1> static constexpr Identidad() requires (r1 == 2u) { return Identidad(std::make_index_sequence<d * d>{}); } TTensor<d, r1> constexpr T() const requires (r1 == 2u) { return T(std::make_index_sequence<d>{}); } // -------------------------------------------------------------------------------------- Operadores // ------------------------------------ Suma (T1 + T2) TTensor<d, r1> constexpr operator +(TTensor<d, r1> const &Tensor) const { return Suma(std::make_index_sequence<d>{}, Tensor); } // ------------------------ Suma-asignación (T1 += T2) TTensor<d, r1> constexpr &operator +=(TTensor<d, r1> const &Tensor) { SumaAsigna(std::make_index_sequence<d>{}, Tensor); return *this; } // ----------------------------------- Resta (T1 - T2) TTensor<d, r1> constexpr operator -(TTensor<d, r1> const &Tensor) const { return Resta(std::make_index_sequence<d>{}, Tensor); } // ----------------------- Resta-asignación (T1 -= T2) TTensor<d, r1> constexpr &operator -=(TTensor<d, r1> const &Tensor) { RestaAsigna(std::make_index_sequence<d>{}, Tensor); return *this; } // ------------------------------------- Negativo (-T) TTensor<d, r1> constexpr operator -() const { return Negativo(std::make_index_sequence<d>{}); } // ----------------------- Producto exterior (T1 * T2) template<std::size_t r2> TTensor<d, r1 + r2> constexpr operator *(TTensor<d, r2> const &Tensor) const { return ProdExt(std::make_index_sequence<d>{}, Tensor); } // -------------- Producto-asignación escalar (T *= s) TTensor<d, r1> constexpr &operator *=(double const Escalar) { ProdAsigna(std::make_index_sequence<d>{}, Escalar); return *this; } // -------------- División-asignación escalar (T /= s) TTensor<d, r1> constexpr &operator /=(double const Escalar) { DivAsigna(std::make_index_sequence<d>{}, Escalar); return *this; } // ----------------------- Producto interior (T1 & T2) template<std::size_t r2> TTensor<d, r1 + r2 - 2u> constexpr operator &(TTensor<d, r2> const &Tensor) const { return ProdInt(std::make_index_sequence<d>{}, Tensor); } // ---------------- Producto interior doble (T1 && T2) template<std::size_t r2> TTensor<d, r1 + r2 - 4u> operator &&(TTensor<d, r2> const &Tensor) const { return ProdIntDbl(std::make_index_sequence<d * d>{}, Tensor); } TTensor<d, r1 - 2u> constexpr operator &&(TTensor<d, 2u> const &Tensor) const requires (r1 == 2u) { return ProdIntDbl(std::make_index_sequence<d>{}, Tensor); } // -------------- Acceso a componentes (T[i][j][k]...) TTensor<d, r1 - 1u> const constexpr &operator [](std::size_t const i) const { return TensorArr[i]; } TTensor<d, r1 - 1u> constexpr &operator [](std::size_t const i) { return TensorArr[i]; } // ------------------------------------------------------------------------------------------ Amigos template<std::size_t, std::size_t> friend struct TTensor; std::istream friend &operator >>(std::istream &is, TTensor<d, r1> &Tensor) { return Tensor.ReadASCII(std::make_index_sequence<d>{}, is); } std::ostream friend &operator <<(std::ostream &os, TTensor<d, r1> const &Tensor) { return Tensor.WriteASCII(std::make_index_sequence<d - 1u>{}, os); } }; // ================================================================================================= // ================================================================================== TTensor<d, 0u> template<std::size_t d> struct TTensor<d, 0u> { double Data; // --------------------------------------------------------------------------------------- Funciones private: template<std::size_t r2, std::size_t... i> TTensor<d, r2> constexpr ProdExt(std::index_sequence<i...>, TTensor<d, r2> const &Tensor) const { return {(Data * Tensor[i])...}; } public: constexpr TTensor() = default; constexpr TTensor(double const Escalar) : Data{Escalar} {} // -------------------------------------------------------------------------------------- Operadores template<std::size_t r2> TTensor<d, r2> constexpr operator *(TTensor<d, r2> const &Tensor) const { return ProdExt(std::make_index_sequence<d>{}, Tensor); } TTensor<d, 0u> constexpr operator *(TTensor<d, 0u> const &Tensor) const { return Data * Tensor.Data; } constexpr operator double() const { return Data; } constexpr operator double &() { return Data; } }; // =========================================================================================== ALIAS // ====================================================================================== TVector<d> template<std::size_t d> using TVector = TTensor<d, 1u>; // ================================================================= OPERADORES Y FUNCIONES EN LÍNEA // =================================================================================== TTensor<d, r> template<std::size_t r, std::size_t d> TTensor<d, r> inline constexpr operator *(TTensor<d, r> const &Tensor, double const Escalar) { return Tensor * TTensor<d, 0u>{Escalar}; } template<std::size_t r, std::size_t d> TTensor<d, r> inline constexpr operator *(double const Escalar, TTensor<d, r> const &Tensor) { return Tensor * TTensor<d, 0u>{Escalar}; } // ================================================================================================= template<std::size_t r, std::size_t d> TTensor<d, r> inline constexpr operator /(TTensor<d, r> const &Tensor, double const Escalar) { return Tensor * TTensor<d, 0u>{1.0 / Escalar}; } // ================================================================================================= template<std::size_t d> double inline mag(TTensor<d, 0u> const &Tensor) { return std::abs(Tensor); } template<std::size_t d> double inline mag(TTensor<d, 1u> const &Tensor) { return std::sqrt(Tensor & Tensor); } template<std::size_t d> double inline mag(TTensor<d, 2u> const &Tensor) { return std::sqrt(Tensor && Tensor); } } // VF namespace VF { // ==================================================================== DECLARACIONES POR ADELANTADO // ========================================================================================== TCelda template<std::size_t> class TCelda; // ================================================================================================= // ======================================================================================= TExprBase template<typename> class TExprBase; // ============================================================================== VARIABLES EN LÍNEA // =========================================================================================== MaxID inline constexpr std::size_t MaxID = std::numeric_limits<std::size_t>::max(); // =========================================================================== DECLARACIÓN DE CLASES // ========================================================================================== TPunto template<std::size_t d> class TPunto final : public TVector<d> { private: std::size_t const ID; // Identificador // --------------------------------------------------------------------------------------- Funciones public: TPunto() = delete; TPunto(std::size_t const ID_) : ID(ID_) {} std::size_t ObtID() const { return ID; } // -------------------------------------------------------------------------------------- Operadores TPunto &operator =(TVector<d> const &Vector) { TVector<d>::operator =(Vector); return *this; } }; // ================================================================================================= // ======================================================================================= TPoligono template<std::size_t d> class TPoligono { private: std::vector<TPunto<d> const *> const PtoPtrVec; // Punteros a puntos en orden antihorario // --------------------------------------------------------------------------------------- Funciones private: template<std::size_t... i> bool IgualFwd(std::index_sequence<i...>, std::size_t const j, TPoligono const &Poligono) const { return (... && (&ObtPunto(i + j) == &Poligono.ObtPunto(i))); } template<std::size_t... i> bool IgualBwd(std::index_sequence<i...>, std::size_t const j, TPoligono const &Poligono) const { return (... && (&ObtPunto(ObtNPunto() + 1u - (i + j)) == &Poligono.ObtPunto(i))); } public: TPoligono() = delete; TPoligono(std::vector<TPunto<d> const *> &&PtoPtrVec_) : PtoPtrVec(std::move(PtoPtrVec_)) {} TPoligono(std::same_as<TPunto<d>> auto const * const... PtoPtr) : PtoPtrVec{PtoPtr...} {} std::size_t ObtNPunto() const { return PtoPtrVec.size(); } TPunto<d> const &ObtPunto(std::size_t const i) const { return *PtoPtrVec[i % ObtNPunto()]; } // -------------------------------------------------------------------------------------- Operadores bool operator ==(TPoligono const &) const; }; // ================================================================================================= // =========================================================================================== TCara template<std::size_t d> class TCara { private: TCelda<d> const * const CeldaPPtr; // Puntero a la celda propietaria TCelda<d> const *CeldaFPtr = nullptr; // Puntero a la celda vecina TPoligono<d> const Poligono; // Puntos que definen a la cara TVector<d> Sf, // Área *dividida* por el volumen de la celda Pf, // Centroide L; // Distancia entre celdas o cara-celda si es CC std::size_t CCID = MaxID; // Identificador del contorno al que pertenece // --------------------------------------------------------------------------------------- Funciones public: TCara() = delete; TCara(TCelda<d> const * const CeldaPPtr_, TPoligono<d> &&Poligono_) : CeldaPPtr(CeldaPPtr_), Poligono(std::move(Poligono_)) {} TCelda<d> const &ObtCeldaP() const { return *CeldaPPtr; } void DefCeldaF(TCelda<d> const * const CeldaFPtr_) { CeldaFPtr = CeldaFPtr_; } TCelda<d> const &ObtCeldaF() const { return *CeldaFPtr; } std::size_t ObtNPunto() const { return Poligono.ObtNPunto(); } TPunto<d> const &ObtPunto(std::size_t const i) const { return Poligono.ObtPunto(i); } void DefCCID(std::size_t const CCID_) { CCID = CCID_; } std::size_t ObtCCID() const { return CCID; } bool EsCC() const { return CCID != MaxID; } void DefSf(); TVector<d> const &ObtSf() const { return Sf; } TVector<d> Calcnf() const { return Sf / mag(Sf); } void DefPf(); TVector<d> const &ObtPf() const { return Pf; } void DefL(); TVector<d> const &ObtL() const { return L; } template<std::size_t r> TTensor<d, r> Interpola(TTensor<d, r> const &φP, TTensor<d, r> const &φF) const { return lerp(φP, φF, ((Pf - ObtCeldaP().ObtP()) & Sf) / (L & Sf)); } // -------------------------------------------------------------------------------------- Operadores bool operator ==(TPoligono<d> const &Poligono) const { return TCara::Poligono == Poligono; } bool operator ==(TCara const &Cara) const { return Poligono == Cara.Poligono; } }; // ================================================================================================= // ========================================================================================== TCelda template<std::size_t d> class TCelda { public: enum ETipo { TRIAN, CUADR, TETRA, HEXA, CUNYA, PIRAM }; // ------------------------------------------------------------------------------------------- Datos private: ETipo const Tipo; // Tipo de celda std::size_t const ID; // Identificador std::vector<TPunto<d> const *> const PtoPtrVec; // Punteros a puntos en orden MSH/VTK std::vector<TCara<d>> CaraVec; // Caras TVector<d> P; // Centroide de la celda double V; // Volumen de la celda std::size_t const GrpID; // Identificador del grupo al que pertenece // --------------------------------------------------------------------------------------- Funciones private: void DefNCara(std::size_t const NCara) { CaraVec.reserve(NCara); } void DefCara(TPoligono<d> &&Poligono) { CaraVec.emplace_back(this, std::move(Poligono)); } public: TCelda() = delete; TCelda(ETipo const, std::size_t const, std::vector<TPunto<d> const *> &&, std::size_t const); ETipo ObtTipo() const { return Tipo; } std::size_t ObtID() const { return ID; } std::size_t ObtNPunto() const { return PtoPtrVec.size(); } TPunto<d> const &ObtPunto(std::size_t const i) const { return *PtoPtrVec[i]; } std::size_t ObtNCara() const { return CaraVec.size(); } TCara<d> const &ObtCara(std::size_t const i) const { return CaraVec[i]; } double ObtV() const { return V; } TVector<d> const &ObtP() const { return P; } std::size_t ObtGrpID() const { return GrpID; } void Inicia(); // ------------------------------------------------------------------------------------------ Amigos auto friend begin(TCelda const &Celda) { return std::begin(Celda.CaraVec); } auto friend begin(TCelda &Celda) { return std::begin(Celda.CaraVec); } auto friend end(TCelda const &Celda) { return std::end(Celda.CaraVec); } auto friend end(TCelda &Celda) { return std::end(Celda.CaraVec); } }; // ================================================================================================= // ========================================================================================== TMalla template<std::size_t d> class TMalla { private: static inline std::vector<TPunto<d>> PuntoVec; // Puntos static inline std::vector<TCelda<d>> CeldaVec; // Celdas static inline std::vector<std::vector<std::size_t>> IDVecVec; // ID de las celdas vecinas static inline std::vector<std::string> GrpStrVec, // Grupos de celdas CCStrVec; // Grupos de caras del contorno // --------------------------------------------------------------------------------------- Funciones public: TMalla() = delete; void static ReadMSH(std::string const &); void static WriteVTK(std::ostream &); std::size_t static ObtNCelda() { return CeldaVec.size(); } TCelda<d> const static &ObtCelda(std::size_t const i) { return CeldaVec[i]; } std::vector<std::size_t> const static &ObtIDVec(std::size_t const i) { return IDVecVec[i]; } std::size_t static ObtNCC() { return CCStrVec.size(); } std::size_t static ObtCCID(std::string_view const); std::size_t static ObtGrpID(std::string_view const); auto static Begin() { return std::cbegin(CeldaVec); } auto static End() { return std::cend(CeldaVec); } }; // =========================================================================================== ALIAS // ======================================================================================== TMalla2D using TMalla2D = TMalla<2u>; // ================================================================================================= // ======================================================================================== TMalla3D using TMalla3D = TMalla<3u>; } // VF namespace VF { // =========================================================================== DECLARACIÓN DE CLASES // =========================================================================================== TCoef template<std::size_t d, std::size_t r> struct TCoef { double aP; std::valarray<double> aF; TTensor<d, r> b; // -------------------------------------------------------------------------------------- Operadores TCoef operator +(TCoef const &Coef) const { return {aP + Coef.aP, aF + Coef.aF, b + Coef.b}; } TCoef operator -(TCoef const &Coef) const { return {aP - Coef.aP, aF - Coef.aF, b - Coef.b}; } TCoef operator -() const { return {-aP, -aF, -b}; } TCoef operator *(TCoef const &Coef) const { return {aP * Coef.aP, aF * Coef.aF, b * Coef.b}; } TCoef operator *(double const Escalar) const { return {aP * Escalar, aF * Escalar, b * Escalar}; } TCoef operator /(double const Escalar) const { return {aP / Escalar, aF / Escalar, b / Escalar}; } }; // ============================================================================= OPERADORES EN LÍNEA // ==================================================================================== TCoef<d, r> template<std::size_t d, std::size_t r> TCoef<d, r> inline operator *(double const Escalar, TCoef<d, r> const &Coef) { return Coef * Escalar; } } // VF namespace VF { // ==================================================================== DECLARACIONES POR ADELANTADO // ======================================================================================= TExprBase template<typename> class TExprBase; // ================================================================================================= // ===================================================================================== TExprUnaria template<std::size_t, std::size_t, typename, typename> class TExprUnaria; // ================================================================================================= // ==================================================================================== TExprBinaria template<std::size_t, std::size_t, typename, typename, typename> class TExprBinaria; // ================================================================================================= // ========================================================================================== TCampo template<std::size_t, std::size_t> class TCampo; // ================================================================================================= // ======================================================================================== TSistema template<std::size_t, std::size_t> class TSistema; // ============================================================================== VARIABLES EN LÍNEA // ========================================================================================== EsExpr template<typename T> inline constexpr bool EsExpr = std::is_base_of_v<TExprBase<T>, T>; template<typename T> inline constexpr bool EsExpr<TExprBase<T>> = true; // ================================================================================================= // ========================================================================================= EsCampo template<typename T> inline constexpr bool EsCampo = false; template<std::size_t d, std::size_t r> inline constexpr bool EsCampo<TCampo<d, r>> = true; // =========================================================================== DECLARACIÓN DE CLASES // ======================================================================================= TExprBase template<typename T> class TExprBase {}; template<std::size_t d, std::size_t r1, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> class TExprBase<T<d, r1, TRest...>> { private: T<d, r1, TRest...> const &Obt() const { return *static_cast<T<d, r1, TRest...> const *>(this); } public: decltype(auto) CalcVal(TCelda<d> const &Celda) const { return Obt().CalcVal(Celda); } decltype(auto) CalcVal(TCara<d> const &Cara) const { return Obt().CalcVal(Cara); } decltype(auto) CalcCoef(TCelda<d> const &Celda) const { return Obt().CalcCoef(Celda); } // -------------------------------------------------------------------------------------- Operadores template<typename... URest, template<std::size_t, std::size_t, typename...> typename U> TExprBinaria<d, r1, T<d, r1, TRest...>, U<d, r1, URest...>, std::plus<>> operator +(TExprBase<U<d, r1, URest...>> const &rhs) const { return {*this, rhs}; } TExprBinaria<d, r1, T<d, r1, TRest...>, TTensor<d, r1>, std::plus<>> operator +(TTensor<d, r1> const &rhs) const { return {*this, rhs}; } template<typename... URest, template<std::size_t, std::size_t, typename...> typename U> TExprBinaria<d, r1, T<d, r1, TRest...>, U<d, r1, URest...>, std::minus<>> operator -(TExprBase<U<d, r1, URest...>> const &rhs) const { return {*this, rhs}; } TExprBinaria<d, r1, T<d, r1, TRest...>, TTensor<d, r1>, std::minus<>> operator -(TTensor<d, r1> const &rhs) const { return {*this, rhs}; } TExprUnaria<d, r1, T<d, r1, TRest...>, std::negate<>> operator -() const { return {*this}; } template<std::size_t r2, typename... URest, template<std::size_t, std::size_t, typename...> typename U> TExprBinaria<d, r1 + r2, T<d, r1, TRest...>, U<d, r2, URest...>, std::multiplies<>> operator *(TExprBase<U<d, r2, URest...>> const &rhs) const { return {*this, rhs}; } template<std::size_t r2> TExprBinaria<d, r1 + r2, T<d, r1, TRest...>, TTensor<d, r2>, std::multiplies<>> operator *(TTensor<d, r2> const &rhs) const { return {*this, rhs}; } TExprBinaria<d, r1, T<d, r1, TRest...>, double, std::multiplies<>> operator *(double const rhs) const { return {*this, rhs}; } template<std::size_t r2, typename... URest, template<std::size_t, std::size_t, typename...> typename U> TExprBinaria<d, r1 + r2 - 2u, T<d, r1, TRest...>, U<d, r2, URest...>, std::bit_and<>> operator &(TExprBase<U<d, r2, URest...>> const &rhs) const { return {*this, rhs}; } template<std::size_t r2> TExprBinaria<d, r1 + r2 - 2u, T<d, r1, TRest...>, TTensor<d, r2>, std::bit_and<>> operator &(TTensor<d, r2> const &rhs) const { return {*this, rhs}; } template<std::size_t r2, typename... URest, template<std::size_t, std::size_t, typename...> typename U> TExprBinaria<d, r1 + r2 - 4u, T<d, r1, TRest...>, U<d, r2, URest...>, std::logical_and<>> operator &&(TExprBase<U<d, r2, URest...>> const &rhs) const { return {*this, rhs}; } template<std::size_t r2> TExprBinaria<d, r1 + r2 - 4u, T<d, r1, TRest...>, TTensor<d, r2>, std::logical_and<>> operator &&(TTensor<d, r2> const &rhs) const { return {*this, rhs}; } template<typename... URest, template<std::size_t, std::size_t, typename...> typename U> TExprBinaria<d, r1, T<d, r1, TRest...>, U<d, 0u, URest...>, std::divides<>> operator /(TExprBase<U<d, 0u, URest...>> const &rhs) const { return {*this, rhs}; } TExprBinaria<d, r1, T<d, r1, TRest...>, double, std::divides<>> operator /(double const rhs) const { return {*this, rhs}; } template<typename... URest, template<std::size_t, std::size_t, typename...> typename U> TSistema<d, r1> operator ==(TExprBase<U<d, r1, URest...>> const &rhs) const { return {*this, rhs}; } TSistema<d, r1> operator ==(TTensor<d, r1> const &rhs) const { return {*this, rhs}; } template<std::false_type = {}> decltype(auto) operator [](std::size_t const i) const requires requires { Obt().template operator []<std::true_type{}>(i); } { return Obt()[i]; } template<std::false_type = {}> decltype(auto) operator [](std::size_t const i) const { return CalcVal(TMalla<d>::ObtCelda(i)); } operator T<d, r1, TRest...> const &() const { return Obt(); } }; // ================================================================================================= // ===================================================================================== TExprUnaria template<std::size_t d, std::size_t r, typename T, typename TOp> class TExprUnaria : public TExprBase<TExprUnaria<d, r, T, TOp>> { public: static constexpr bool EsIndexada = true; private: std::conditional_t<EsCampo<T>, T const &, T const> rhs; [[no_unique_address]] TOp const Op = {}; // --------------------------------------------------------------------------------------- Funciones public: TExprUnaria(T const &rhs_) : rhs(rhs_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return Op(rhs.CalcVal(Celda)); } TTensor<d, r> CalcVal(TCara<d> const &Cara) const { return Op(rhs.CalcVal(Cara)); } decltype(auto) CalcCoef(TCelda<d> const &Celda) const { return Op(rhs.CalcCoef(Celda)); } // -------------------------------------------------------------------------------------- Operadores template<std::true_type = {}> TTensor<d, r> operator [](std::size_t const i) const { return Op(rhs[i]); } }; // ================================================================================================= // ==================================================================================== TExprBinaria template<std::size_t d, std::size_t r, typename T, typename U, typename TOp> class TExprBinaria : public TExprBase<TExprBinaria<d, r, T, U, TOp>> { private: std::conditional_t<EsCampo<T>, T const &, T const> lhs; std::conditional_t<EsCampo<U>, U const &, U const> rhs; [[no_unique_address]] TOp const Op = {}; // --------------------------------------------------------------------------------------- Funciones public: TExprBinaria(T const &lhs_, U const &rhs_) : lhs(lhs_), rhs(rhs_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const requires (EsExpr<T> && EsExpr<U>) { return Op(lhs.CalcVal(Celda), rhs.CalcVal(Celda)); } TTensor<d, r> CalcVal(TCelda<d> const &Celda) const requires (EsExpr<T> && !EsExpr<U>) { return Op(lhs.CalcVal(Celda), rhs); } TTensor<d, r> CalcVal(TCelda<d> const &Celda) const requires (!EsExpr<T> && EsExpr<U>) { return Op(lhs, rhs.CalcVal(Celda)); } TTensor<d, r> CalcVal(TCara<d> const &Cara) const requires (EsExpr<T> && EsExpr<U>) { return Op(lhs.CalcVal(Cara), rhs.CalcVal(Cara)); } TTensor<d, r> CalcVal(TCara<d> const &Cara) const requires (EsExpr<T> && !EsExpr<U>) { return Op(lhs.CalcVal(Cara), rhs); } TTensor<d, r> CalcVal(TCara<d> const &Cara) const requires (!EsExpr<T> && EsExpr<U>) { return Op(lhs, rhs.CalcVal(Cara)); } decltype(auto) CalcCoef(TCelda<d> const &Celda) const requires (EsExpr<T> && EsExpr<U>) { return Op(lhs.CalcCoef(Celda), rhs.CalcCoef(Celda)); } decltype(auto) CalcCoef(TCelda<d> const &Celda) const requires (EsExpr<T> && !EsExpr<U>) { return Op(lhs.CalcCoef(Celda), rhs); } decltype(auto) CalcCoef(TCelda<d> const &Celda) const requires (!EsExpr<T> && EsExpr<U>) { return Op(lhs, rhs.CalcCoef(Celda)); } // -------------------------------------------------------------------------------------- Operadores template<std::true_type = std::true_type{}> TTensor<d, r> operator [](std::size_t const i) const requires (EsExpr<T> && EsExpr<U>) { return Op(lhs[i], rhs[i]); } template<std::true_type = std::true_type{}> TTensor<d, r> operator [](std::size_t const i) const requires (EsExpr<T> && !EsExpr<U>) { return Op(lhs[i], rhs); } template<std::true_type = std::true_type{}> TTensor<d, r> operator [](std::size_t const i) const requires (!EsExpr<T> && EsExpr<U>) { return Op(lhs, rhs[i]); } // ------------------------------------------------------------------------------------------ Amigos template<std::size_t, std::size_t, typename> friend class TExprDiv; template<std::size_t, std::size_t, typename> friend class TExprLap; template<std::size_t, std::size_t, typename> friend class TExprSp; }; // ============================================================================= OPERADORES EN LÍNEA // ======================================================================================= TExprBase template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, r, TTensor<d, r>, T<d, r, TRest...>, std::plus<>> inline operator +(TTensor<d, r> const &lhs, TExprBase<T<d, r, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, 0u, T<d, 0u, TRest...>, double, std::plus<>> inline operator +(TExprBase<T<d, 0u, TRest...>> const &lhs, double const rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, 0u, double, T<d, 0u, TRest...>, std::plus<>> inline operator +(double const lhs, TExprBase<T<d, 0u, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, r, TTensor<d, r>, T<d, r, TRest...>, std::minus<>> inline operator -(TTensor<d, r> const &lhs, TExprBase<T<d, r, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, 0u, T<d, 0u, TRest...>, double, std::minus<>> inline operator -(TExprBase<T<d, 0u, TRest...>> const &lhs, double const rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, 0u, double, T<d, 0u, TRest...>, std::minus<>> inline operator -(double const lhs, TExprBase<T<d, 0u, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, std::size_t r1, std::size_t r2, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, r1 + r2, TTensor<d, r1>, T<d, r2, TRest...>, std::multiplies<>> inline operator *(TTensor<d, r1> const &lhs, TExprBase<T<d, r2, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, r, double, T<d, r, TRest...>, std::multiplies<>> inline operator *(double const lhs, TExprBase<T<d, r, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, std::size_t r1, std::size_t r2, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, r1 + r2 - 2u, TTensor<d, r1>, T<d, r2, TRest...>, std::bit_and<>> inline operator &(TTensor<d, r1> const &lhs, TExprBase<T<d, r2, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, std::size_t r1, std::size_t r2, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, r1 + r2 - 4u, TTensor<d, r1>, T<d, r2, TRest...>, std::logical_and<>> inline operator &&(TTensor<d, r1> const &lhs, TExprBase<T<d, r2, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, r, TTensor<d, r>, T<d, 0u, TRest...>, std::divides<>> inline operator /(TTensor<d, r> const &lhs, TExprBase<T<d, 0u, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<std::size_t d, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<d, 0u, double, T<d, 0u, TRest...>, std::divides<>> inline operator /(double const lhs, TExprBase<T<d, 0u, TRest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T, typename... URest, template<std::size_t, std::size_t, typename...> typename U> TExprBinaria<3u, 1u, T<3u, 1u, TRest...>, U<3u, 1u, URest...>, std::bit_xor<>> inline operator ^(TExprBase<T<3u, 1u, TRest...>> const &lhs, TExprBase<U<3u, 1u, URest...>> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<3u, 1u, T<3u, 1u, TRest...>, TTensor<3u, 1u>, std::bit_xor<>> inline operator ^(TExprBase<T<3u, 1u, TRest...>> const &lhs, TTensor<3u, 1u> const &rhs) { return {lhs, rhs}; } // ================================================================================================= template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprBinaria<3u, 1u, TTensor<3u, 1u>, T<3u, 1u, TRest...>, std::bit_xor<>> inline operator ^(TTensor<3u, 1u> const &lhs, TExprBase<T<3u, 1u, TRest...>> const &rhs) { return {lhs, rhs}; } } // VF namespace VF { // ================================================================================================= // ============================================================================================= mag template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T, typename TOp = decltype([](TTensor<d, r> const &Tensor) { return mag(Tensor); })> TExprUnaria<d, 0u, T<d, r, TRest...>, TOp> inline mag(TExprBase<T<d, r, TRest...>> const &Expr) { return {Expr}; } // ================================================================================================= // ============================================================================================= sum template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TTensor<d, r> inline sum(TExprBase<T<d, r, TRest...>> const &Expr) { TTensor<d, r> Sum = {}; #pragma omp declare reduction(Op: TTensor<d, r>: omp_out += omp_in) initializer(omp_priv = {}) #pragma omp parallel for reduction(Op: Sum) for (std::size_t i = 0u; i < TMalla<d>::ObtNCelda(); ++i) Sum += Expr[i]; return Sum; } } // VF namespace VF { // =========================================================================== DECLARACIÓN DE CLASES // ================================================================================== TEsCCExplicita template<typename T> struct TEsCCExplicita : std::false_type {}; template<typename T> inline constexpr bool EsCCExplicita = TEsCCExplicita<T>::value; // ================================================================================================= // ========================================================================================= TCCBase template<std::size_t d, std::size_t r> class TCCBase { public: virtual ~TCCBase() = default; std::tuple<double, TTensor<d, r>> virtual CalcVal(TCara<d> const &) const; std::tuple<double, TTensor<d, r>> virtual CalcGrad(TCara<d> const &) const; }; // ================================================================================================= // ====================================================================================== TDirichlet template<std::size_t d, std::size_t r> class TDirichlet : public TCCBase<d, r> { private: TTensor<d, r> const φb = {}; // φf = φb // --------------------------------------------------------------------------------------- Funciones public: TDirichlet() = default; TDirichlet(TTensor<d, r> const &φb_) : φb(φb_) {} std::tuple<double, TTensor<d, r>> virtual CalcVal(TCara<d> const &) const override { return {0.0, φb}; } std::tuple<double, TTensor<d, r>> virtual CalcGrad(TCara<d> const &Cara) const override { double const magL = mag(Cara.ObtL()); return {-1.0 / magL, φb / magL}; } }; // ================================================================================================= // ======================================================================================== TNeumann template<std::size_t d, std::size_t r> class TNeumann : public TCCBase<d, r> { private: TTensor<d, r> const gb = {}; // nf & ∇ * φf = gb // --------------------------------------------------------------------------------------- Funciones public: TNeumann() = default; TNeumann(TTensor<d, r> const &gb_) : gb(gb_) {} std::tuple<double, TTensor<d, r>> virtual CalcVal(TCara<d> const &Cara) const override { return {1.0, mag(Cara.ObtL())*gb}; } std::tuple<double, TTensor<d, r>> virtual CalcGrad(TCara<d> const &) const override { return {0.0, gb}; } }; // ================================================================================================= // ========================================================================================== TRobin template<std::size_t d, std::size_t r> class TRobin : public TCCBase<d, r> { private: double const k, h; TTensor<d, r> const hφb; // k * (nf & ∇ * φf) + h * φf = h * φb // --------------------------------------------------------------------------------------- Funciones public: TRobin() = delete; TRobin(double const k_, double const h_, TTensor<d, r> const &φb) : k(k_), h(h_), hφb(h_ * φb) {} std::tuple<double, TTensor<d, r>> virtual CalcGrad(TCara<d> const &Cara) const override { return {-h / (k + h * mag(Cara.ObtL())), (1.0 / (k + h * mag(Cara.ObtL()))) * hφb}; } }; // ======================================================================== IMPLEMENTACIÓN DE CLASES // ========================================================================================= TCCBase template<std::size_t d, std::size_t r> std::tuple<double, TTensor<d, r>> TCCBase<d, r>::CalcVal(TCara<d> const &Cara) const { auto const [aP, b] = CalcGrad(Cara); double const magL = mag(Cara.ObtL()); return {1.0 + magL * aP, magL * b}; } // ================================================================================================= template<std::size_t d, std::size_t r> std::tuple<double, TTensor<d, r>> TCCBase<d, r>::CalcGrad(TCara<d> const &Cara) const { auto const [aP, b] = CalcVal(Cara); double const magL = mag(Cara.ObtL()); return {(aP - 1.0) / magL, b / magL}; } } // VF namespace VF { // =========================================================================== DECLARACIÓN DE CLASES // ========================================================================================== TCampo template<std::size_t d, std::size_t r> class TCampo : public TExprBase<TCampo<d, r>> { private: class TCCPtr { private: std::unique_ptr<TCCBase<d, r>> CCPtr = std::make_unique<TNeumann<d, r>>(); public: TCCPtr() = default; // Por defecto Neumann homogénea TCCPtr &operator =(std::unique_ptr<TCCBase<d, r>> &&CCPtr_) { CCPtr = std::move(CCPtr_); return *this; } operator TCCBase<d, r> const *() const { return CCPtr.get(); } }; // ------------------------------------------------------------------------------------------- Datos private: std::size_t const NCelda = TMalla<d>::ObtNCelda(); TTensor<d, r> *TensorPtr = new TTensor<d, r>[NCelda]{}; std::vector<TCCPtr> CCPtrVec = std::vector<TCCPtr>(TMalla<d>::ObtNCC()); // --------------------------------------------------------------------------------------- Funciones private: void Copia(TCampo const &Campo) { std::copy_n(Campo.TensorPtr, NCelda, TensorPtr); } void Llena(TTensor<d, r> const &Tensor) { std::fill_n(TensorPtr, NCelda, Tensor); } template<typename T> void Asigna(TExprBase<T> const &); template<typename T, typename... TArgs> void DefCC(std::false_type, std::string_view const CCStr, TArgs &&...Args) { CCPtrVec[TMalla<d>::ObtCCID(CCStr)] = std::make_unique<T>(std::forward<TArgs>(Args)...); } template<typename T, typename... TArgs> void DefCC(std::true_type, std::string_view const CCStr, TArgs &&...Args) { DefCC<T>(std::false_type{}, CCStr, *this, std::forward<TArgs>(Args)...); } TCCBase<d, r> const *ObtCCPtr(TCara<d> const &Cara) const { return CCPtrVec[Cara.ObtCCID()]; } public: TCampo() = default; ~TCampo() { delete[] TensorPtr; } TCampo(TCampo const &Campo) { Copia(Campo); } TCampo(TCampo &&Campo) : TensorPtr(std::exchange(Campo.TensorPtr, nullptr)) {} TCampo(TTensor<d, r> const &Tensor) { Llena(Tensor); } template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TCampo(TExprBase<T<d, r, TRest...>> const &Expr) { Asigna(Expr); } template<template<std::size_t, std::size_t> typename T, typename... TArgs> void DefCC(std::string_view const CCStr, TArgs &&...Args) { DefCC<T<d, r>>(TEsCCExplicita<T<d, r>>{}, CCStr, std::forward<TArgs>(Args)...); } template<template<std::size_t, std::size_t, typename...> typename T, typename... TArgs, typename U = T<d, r, std::remove_cvref_t<TArgs>...>> requires (sizeof...(TArgs) != 0u) void DefCC(std::string_view const CCStr, TArgs &&...Args) { DefCC<U>(TEsCCExplicita<U>{}, CCStr, std::forward<TArgs>(Args)...); } TTensor<d, r> const &CalcVal(TCelda<d> const &Celda) const { return TensorPtr[Celda.ObtID()]; } TTensor<d, r> CalcVal(TCara<d> const &) const; std::tuple<double, TTensor<d, r>> CalcValCC(TCara<d> const &Cara) const { return ObtCCPtr(Cara)->CalcVal(Cara); } TTensor<d, r + 1u> CalcGrad(TCelda<d> const &) const; TTensor<d, r> CalcGrad(TCara<d> const &) const; std::tuple<double, TTensor<d, r>> CalcGradCC(TCara<d> const &Cara) const { return ObtCCPtr(Cara)->CalcGrad(Cara); } TTensor<d, r + 1u> CalcGradT(TCelda<d> const &) const; std::conditional_t<r != 0u, TTensor<d, r - 1u>, void> CalcDiv(TCelda<d> const &) const; TTensor<d, r> CalcDiv(TTensor<d, 1u> const &, TCelda<d> const &) const; template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TTensor<d, r> CalcDiv(TExprBase<T<d, 1u, TRest...>> const &, TCelda<d> const &) const; TTensor<d, r> CalcRot(TCelda<d> const &) const = delete; TTensor<d, r> CalcLap(TCelda<d> const &) const; template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TTensor<d, r> CalcLap(TExprBase<T<d, 0u, TRest...>> const &, TCelda<d> const &) const; TTensor<d, r> CalcLapT(TCelda<d> const &) const; TTensor<d, r> const &CalcCoef(TCelda<d> const &Celda) const { return TensorPtr[Celda.ObtID()]; } void WriteVTK(std::ostream &, std::string_view const) const; // -------------------------------------------------------------------------------------- Operadores TCampo &operator =(TCampo const &Campo) { Copia(Campo); return *this; } TCampo &operator =(TCampo &&Campo) noexcept { delete[] TensorPtr; TensorPtr = std::exchange(Campo.TensorPtr, nullptr); return *this; } TCampo &operator =(TTensor<d, r> const &Tensor) { Llena(Tensor); return *this; } template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TCampo &operator =(TExprBase<T<d, r, TRest...>> const &Expr) { Asigna(Expr); return *this; } template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TCampo &operator +=(TExprBase<T<d, r, TRest...>> const &); TCampo &operator +=(TTensor<d, r> const &); template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TCampo &operator -=(TExprBase<T<d, r, TRest...>> const &); TCampo &operator -=(TTensor<d, r> const &); template<std::true_type = std::true_type{}> TTensor<d, r> const &operator [](std::size_t const i) const { return TensorPtr[i]; } TTensor<d, r> &operator [](std::size_t const i) { return TensorPtr[i]; } // ------------------------------------------------------------------------------------------ Amigos TTensor<d, r> const friend *begin(TCampo const &Campo) { return Campo.TensorPtr; } TTensor<d, r> friend *begin(TCampo &Campo) { return Campo.TensorPtr; } TTensor<d, r> const friend *end(TCampo const &Campo) { return Campo.TensorPtr + Campo.NCelda; } TTensor<d, r> friend *end(TCampo &Campo) { return Campo.TensorPtr + Campo.NCelda; } }; // ================================================================================================= template<> void TCampo<2u, 0u>::CalcDiv(TCelda<2u> const &) const = delete; // ================================================================================================= template<> void TCampo<3u, 0u>::CalcDiv(TCelda<3u> const &) const = delete; // ================================================================================================= template<> TTensor<2u, 0u> TCampo<2u, 0u>::CalcLapT(TCelda<2u> const &) const = delete; // ================================================================================================= template<> TTensor<3u, 0u> TCampo<3u, 0u>::CalcLapT(TCelda<3u> const &) const = delete; // ================================================================================================= template<> TTensor<3u, 1u> TCampo<3u, 1u>::CalcRot(TCelda<3u> const &) const; // =========================================================================================== ALIAS // =================================================================================== TCampoEscalar template<std::size_t d> using TCampoEscalar = TCampo<d, 0u>; // ================================================================================================= // ================================================================================= TCampoEscalar2D using TCampoEscalar2D = TCampoEscalar<2u>; // ================================================================================================= // ================================================================================= TCampoEscalar3D using TCampoEscalar3D = TCampoEscalar<3u>; // ================================================================================================= // ================================================================================= TCampoVectorial template<std::size_t d> using TCampoVectorial = TCampo<d, 1u>; // ================================================================================================= // =============================================================================== TCampoVectorial2D using TCampoVectorial2D = TCampoVectorial<2u>; // ================================================================================================= // =============================================================================== TCampoVectorial3D using TCampoVectorial3D = TCampoVectorial<3u>; // ======================================================================== IMPLEMENTACIÓN DE CLASES // ========================================================================================== TCampo template<std::size_t d, std::size_t r> template<typename T> void TCampo<d, r>::Asigna(TExprBase<T> const &Expr) { #pragma omp parallel for for (std::size_t i = 0u; i < NCelda; ++i) TensorPtr[i] = Expr[i]; } // ================================================================================================= template<std::size_t d, std::size_t r> TCampo<d, r> &TCampo<d, r>::operator +=(TTensor<d, r> const &Tensor) { #pragma omp parallel for for (std::size_t i = 0u; i < NCelda; ++i) TensorPtr[i] += Tensor; return *this; } // ================================================================================================= template<std::size_t d, std::size_t r> template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TCampo<d, r> &TCampo<d, r>::operator +=(TExprBase<T<d, r, TRest...>> const &Expr) { #pragma omp parallel for for (std::size_t i = 0u; i < NCelda; ++i) TensorPtr[i] += Expr[i]; return *this; } // ================================================================================================= template<std::size_t d, std::size_t r> TCampo<d, r> &TCampo<d, r>::operator -=(TTensor<d, r> const &Tensor) { #pragma omp parallel for for (std::size_t i = 0u; i < NCelda; ++i) TensorPtr[i] -= Tensor; return *this; } // ================================================================================================= template<std::size_t d, std::size_t r> template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TCampo<d, r> &TCampo<d, r>::operator -=(TExprBase<T<d, r, TRest...>> const &Expr) { #pragma omp parallel for for (std::size_t i = 0u; i < NCelda; ++i) TensorPtr[i] -= Expr[i]; return *this; } // ================================================================================================= template<std::size_t d, std::size_t r> TTensor<d, r> TCampo<d, r>::CalcVal(TCara<d> const &Cara) const { if (Cara.EsCC()) [[unlikely]] { auto const [aP, b] = CalcValCC(Cara); return aP * CalcVal(Cara.ObtCeldaP()) + b; } return Cara.Interpola(CalcVal(Cara.ObtCeldaP()), CalcVal(Cara.ObtCeldaF())); } // ================================================================================================= template<std::size_t d, std::size_t r> TTensor<d, r + 1u> TCampo<d, r>::CalcGrad(TCelda<d> const &Celda) const { TTensor<d, r + 1u> Grad = {}; for (auto &Cara : Celda) Grad += Cara.ObtSf() * CalcVal(Cara); return Grad; } // ================================================================================================= template<std::size_t d, std::size_t r> TTensor<d, r> TCampo<d, r>::CalcGrad(TCara<d> const &Cara) const { if (Cara.EsCC()) [[unlikely]] { auto const [aP, b] = CalcGradCC(Cara); return aP * CalcVal(Cara.ObtCeldaP()) + b; } return Cara.Calcnf() & Cara.Interpola(CalcGrad(Cara.ObtCeldaP()), CalcGrad(Cara.ObtCeldaF())); } // ================================================================================================= template<std::size_t d, std::size_t r> TTensor<d, r + 1u> TCampo<d, r>::CalcGradT(TCelda<d> const &Celda) const { TTensor<d, r + 1u> GradT = {}; for (auto &Cara : Celda) GradT += CalcVal(Cara) * Cara.ObtSf(); return GradT; } // ================================================================================================= template<std::size_t d, std::size_t r> std::conditional_t<r != 0u, TTensor<d, r - 1u>, void> TCampo<d, r>::CalcDiv(TCelda<d> const &Celda) const { TTensor<d, r - 1u> Div = {}; for (auto &Cara : Celda) Div += Cara.ObtSf() & CalcVal(Cara); return Div; } // ================================================================================================= template<std::size_t d, std::size_t r> TTensor<d, r> TCampo<d, r>::CalcDiv(TTensor<d, 1u> const &U, TCelda<d> const &Celda) const { TTensor<d, r> Div = {}; for (auto &Cara : Celda) Div += (Cara.ObtSf() & U) * CalcVal(Cara); return Div; } // ================================================================================================= template<std::size_t d, std::size_t r> template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TTensor<d, r> TCampo<d, r>::CalcDiv(TExprBase<T<d, 1u, TRest...>> const &U, TCelda<d> const &Celda) const { TTensor<d, r> Div = {}; for (auto &Cara : Celda) Div += (Cara.ObtSf() & U.CalcVal(Cara)) * CalcVal(Cara); return Div; } // ================================================================================================= template<std::size_t d, std::size_t r> TTensor<d, r> TCampo<d, r>::CalcLap(TCelda<d> const &Celda) const { TTensor<d, r + 1u> const Grad = CalcGrad(Celda); TTensor<d, r> Lap = {}; for (auto &Cara : Celda) if (Cara.EsCC()) [[unlikely]] { auto const [aP, b] = CalcGradCC(Cara); Lap += mag(Cara.ObtSf()) * (aP * CalcVal(Celda) + b); } else Lap += Cara.ObtSf() & Cara.Interpola(Grad, CalcGrad(Cara.ObtCeldaF())); return Lap; } // ================================================================================================= template<std::size_t d, std::size_t r> template<typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TTensor<d, r> TCampo<d, r>::CalcLap(TExprBase<T<d, 0u, TRest...>> const &Γ, TCelda<d> const &Celda) const { TTensor<d, r + 1u> const Grad = CalcGrad(Celda); TTensor<d, r> Lap = {}; for (auto &Cara : Celda) { TTensor<d, 0u> Γf = Γ.CalcVal(Cara); if (Cara.EsCC()) [[unlikely]] { auto const [aP, b] = CalcGradCC(Cara); Lap += Γf * mag(Cara.ObtSf()) * (aP * CalcVal(Cara.ObtCeldaP()) + b); } else Lap += Γf * (Cara.ObtSf() & Cara.Interpola(Grad, CalcGrad(Cara.ObtCeldaF()))); } return Lap; } // ================================================================================================= template<std::size_t d, std::size_t r> TTensor<d, r> TCampo<d, r>::CalcLapT(TCelda<d> const &Celda) const { TTensor<d, r + 1u> const GradT = CalcGradT(Celda); TTensor<d, r> LapT = {}; for (auto &Cara : Celda) if (Cara.EsCC()) [[unlikely]] { TVector<d> const nf = Cara.Calcnf(); auto const [aP, b] = CalcGradCC(Cara); LapT += Cara.ObtSf() & (GradT + (aP * CalcVal(Celda) + b - (GradT & nf)) * nf); } else LapT += Cara.ObtSf() & Cara.Interpola(GradT, CalcGradT(Cara.ObtCeldaF())); return LapT; } // ================================================================================================= template<std::size_t d, std::size_t r> void TCampo<d, r>::WriteVTK(std::ostream &os, std::string_view const Nombre) const { if constexpr (r == 0u) { os << "SCALARS " << Nombre << " double 1" << std::endl; os << "LOOKUP_TABLE default" << std::endl; for (auto &Tensor : *this) os << Tensor << std::endl; } else if constexpr (r == 1u) { os << "VECTORS " << Nombre << " double" << std::endl; for (auto &Tensor : *this) if constexpr (d == 2u) os << Tensor[0u] << ' ' << Tensor[1u] << " 0" << std::endl; else if constexpr (d == 3u) os << Tensor; } else if constexpr (r == 2u) { os << "TENSORS " << Nombre << " double" << std::endl; for (auto &Tensor : *this) if constexpr (d == 2u) { os << Tensor[0u][0u] << ' ' << Tensor[0u][1u] << " 0" << std::endl; os << Tensor[1u][0u] << ' ' << Tensor[1u][1u] << " 0" << std::endl; os << "0 0 0" << std::endl; } else if constexpr (d == 3u) os << Tensor; } } } // VF // FUNCIÓN OPERACIÓN IMPLÍCITO // ------- --------- --------- // ∂φ // d(φ) / dt ---- Sí // ∂t // // div(φ) ∇ & φ No // // div(U * φ) ∇ & (U * φ) Sí // // grad(φ) ∇ * φ No // // T // gradT(φ) (∇ * φ) No // // lap(φ) ∇ & ∇ * φ Sí // // div(Γ * grad(φ)) ∇ & (Γ * ∇ * φ) Sí // // T // lapT(φ) ∇ & (∇ * φ) No // // rot(φ) ∇ ^ φ No // // Sp(φ) φ Sí // // Sp(ρ * φ) ρ * φ Sí namespace VF { // =========================================================================== DECLARACIÓN DE CLASES // ========================================================================================== TExprd template<std::size_t d, std::size_t r> class TExprd : public TExprBase<TExprd<d, r>> { private: TCampo<d, r> const &φ; // --------------------------------------------------------------------------------------- Funciones public: TExprd(TCampo<d, r> const &φ_) : φ(φ_) {} TTensor<d, r> CalcVal(TCelda<d> const &) const = delete; TTensor<d, r> CalcVal(TCara<d> const &) const = delete; TCoef<d, r> CalcCoef(TCelda<d> const &Celda) const { return {1.0, std::valarray<double>(Celda.ObtNCara()), -φ.CalcVal(Celda)}; } }; // ================================================================================================= // ======================================================================================== TExprDiv template<std::size_t d, std::size_t r, typename T = void> class TExprDiv : public TExprBase<TExprDiv<d, r, T>> { private: std::conditional_t<EsCampo<T>, T const &, T const> U; TCampo<d, r> const &φ; // --------------------------------------------------------------------------------------- Funciones private: double CalcSfUf(TCara<d> const &Cara) const requires EsExpr<T> { return Cara.ObtSf() & U.CalcVal(Cara); } double CalcSfUf(TCara<d> const &Cara) const { return Cara.ObtSf() & U; } std::tuple<double, double, TTensor<d, r>> CalcCoef(TCara<d> const &, double const) const; public: TExprDiv(TExprBinaria<d, r + 1u, T, TCampo<d, r>, std::multiplies<>> const &Expr) : U(Expr.lhs), φ(Expr.rhs) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return φ.CalcDiv(U, Celda); } TCoef<d, r> CalcCoef(TCelda<d> const &) const; }; // ================================================================================================= // ============================================================================ TExprDiv<d, r, void> template<std::size_t d, std::size_t r> class TExprDiv<d, r, void> : public TExprBase<TExprDiv<d, r, void>> { private: TCampo<d, r + 1u> const &φ; // --------------------------------------------------------------------------------------- Funciones public: TExprDiv(TCampo<d, r + 1u> const &φ_) : φ(φ_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return φ.CalcDiv(Celda); } TTensor<d, r> CalcCoef(TCelda<d> const &Celda) const { return CalcVal(Celda); } }; // ================================================================================================= // ======================================================================================= TExprGrad template<std::size_t d, std::size_t r> class TExprGrad : public TExprBase<TExprGrad<d, r>> { private: TCampo<d, r - 1u> const &φ; // --------------------------------------------------------------------------------------- Funciones public: TExprGrad(TCampo<d, r - 1u> const &φ_) : φ(φ_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return φ.CalcGrad(Celda); } TTensor<d, r> CalcCoef(TCelda<d> const &Celda) const { return CalcVal(Celda); } // ------------------------------------------------------------------------------------------ Amigos template<std::size_t, std::size_t, typename> friend class TExprLap; }; // ================================================================================================= // ====================================================================================== TExprGradT template<std::size_t d, std::size_t r> class TExprGradT : public TExprBase<TExprGradT<d, r>> { private: TCampo<d, r - 1u> const &φ; // --------------------------------------------------------------------------------------- Funciones public: TExprGradT(TCampo<d, r - 1u> const &φ_) : φ(φ_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return φ.CalcGradT(Celda); } TTensor<d, r> CalcCoef(TCelda<d> const &Celda) const { return CalcVal(Celda); } }; // ================================================================================================= // ==================================================================================== TExprLapBase template<std::size_t d, std::size_t r> class TExprLapBase { protected: TCampo<d, r> const &φ; // --------------------------------------------------------------------------------------- Funciones protected: TExprLapBase(TCampo<d, r> const &φ_) : φ(φ_) {} std::tuple<double, double, TTensor<d, r>> CalcCoef(TCara<d> const &, TTensor<d, r + 1u> const &) const; }; // ================================================================================================= // ======================================================================================== TExprLap template<std::size_t d, std::size_t r, typename T = void> class TExprLap : public TExprLapBase<d, r>, public TExprBase<TExprLap<d, r, T>> { private: using TExprLapBase<d, r>::φ; std::conditional_t<EsCampo<T>, T const &, T const> Γ; // --------------------------------------------------------------------------------------- Funciones public: TExprLap(TExprBinaria<d, r + 1u, T, TExprGrad<d, r + 1u>, std::multiplies<>> const &Expr) : TExprLapBase<d, r>(Expr.rhs.φ), Γ(Expr.lhs) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return φ.CalcLap(Γ, Celda); } TCoef<d, r> CalcCoef(TCelda<d> const &) const; }; // ================================================================================================= // ============================================================================ TExprLap<d, r, void> template<std::size_t d, std::size_t r> class TExprLap<d, r, void> : public TExprLapBase<d, r>, public TExprBase<TExprLap<d, r, void>> { private: using TExprLapBase<d, r>::φ; // --------------------------------------------------------------------------------------- Funciones public: TExprLap(TCampo<d, r> const &φ_) : TExprLapBase<d, r>(φ_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return φ.CalcLap(Celda); } TCoef<d, r> CalcCoef(TCelda<d> const &) const; }; // ================================================================================================= // ======================================================================================= TExprLapT template<std::size_t d, std::size_t r> class TExprLapT : public TExprBase<TExprLapT<d, r>> { private: TCampo<d, r> const &φ; // --------------------------------------------------------------------------------------- Funciones public: TExprLapT(TCampo<d, r> const &φ_) : φ(φ_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return φ.CalcLapT(Celda); } TTensor<d, r> CalcCoef(TCelda<d> const &Celda) const { return CalcVal(Celda); } }; // ================================================================================================= // ======================================================================================== TExprRot template<std::size_t d, std::size_t r> class TExprRot : public TExprBase<TExprRot<d, r>> { private: TCampo<d, r> const &φ; // --------------------------------------------------------------------------------------- Funciones public: TExprRot(TCampo<d, r> const &φ_) : φ(φ_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return φ.CalcRot(Celda); } TTensor<d, r> CalcCoef(TCelda<d> const &Celda) const { return CalcVal(Celda); } }; // ================================================================================================= // ========================================================================================= TExprSp template<std::size_t d, std::size_t r, typename T = void> class TExprSp : public TExprBase<TExprSp<d, r, T>> { private: std::conditional_t<EsCampo<T>, T const &, T const> ρ; TCampo<d, r> const &φ; // --------------------------------------------------------------------------------------- Funciones public: TExprSp(TExprBinaria<d, r, T, TCampo<d, r>, std::multiplies<>> const &Expr) : ρ(Expr.lhs), φ(Expr.rhs) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return ρ.CalcVal(Celda) * φ.CalcVal(Celda); } TCoef<d, r> CalcCoef(TCelda<d> const &) const; }; // ================================================================================================= // ============================================================================ TExprSp<d, r, void> template<std::size_t d, std::size_t r> class TExprSp<d, r, void> : public TExprBase<TExprSp<d, r, void>> { private: TCampo<d, r> const &φ; // --------------------------------------------------------------------------------------- Funciones public: TExprSp(TCampo<d, r> const &φ_) : φ(φ_) {} TTensor<d, r> const &CalcVal(TCelda<d> const &Celda) const { return φ.CalcVal(Celda); } TCoef<d, r> CalcCoef(TCelda<d> const &Celda) const { return {1.0, std::valarray<double>(Celda.ObtNCara()), {}}; } }; // ======================================================================== IMPLEMENTACIÓN DE CLASES // ======================================================================================== TExprDiv template<std::size_t d, std::size_t r, typename T> std::tuple<double, double, TTensor<d, r>> TExprDiv<d, r, T>::CalcCoef(TCara<d> const &Cara, double const SfUf) const { if (SfUf > 0.0) return {SfUf, 0.0, {}}; if (Cara.EsCC()) [[unlikely]] { auto const [aP, b] = φ.CalcValCC(Cara); return {SfUf * aP, 0.0, SfUf * b}; } return {0.0, SfUf, {}}; } // ================================================================================================= template<std::size_t d, std::size_t r, typename T> TCoef<d, r> TExprDiv<d, r, T>::CalcCoef(TCelda<d> const &Celda) const { std::size_t const NCara = Celda.ObtNCara(); double aP = 0.0; std::valarray<double> aF(NCara); TTensor<d, r> b = {}; for (std::size_t i = 0u; i < NCara; ++i) { TCara<d> const &Cara = Celda.ObtCara(i); auto const [aP_, aF_, b_] = CalcCoef(Cara, CalcSfUf(Cara)); aP += aP_; aF[i] = aF_; b += b_; } return {aP, aF, b}; } // ================================================================================================= // ==================================================================================== TExprLapBase template<std::size_t d, std::size_t r> std::tuple<double, double, TTensor<d, r>> TExprLapBase<d, r>::CalcCoef(TCara<d> const &Cara, TTensor<d, r + 1u> const &Grad) const { TVector<d> const &Sf = Cara.ObtSf(); if (Cara.EsCC()) [[unlikely]] { auto const [aP, b] = φ.CalcGradCC(Cara); return {mag(Sf) * aP, 0.0, mag(Sf) * b}; } TVector<d> const &L = Cara.ObtL(), SfPar = ((Sf & Sf) / (L & Sf)) * L; double const aF = mag(SfPar) / mag(L); return {-aF, aF, (Sf - SfPar) & Cara.Interpola(Grad, φ.CalcGrad(Cara.ObtCeldaF()))}; } // ================================================================================================= // ======================================================================================== TExprLap template<std::size_t d, std::size_t r, typename T> TCoef<d, r> TExprLap<d, r, T>::CalcCoef(TCelda<d> const &Celda) const { std::size_t const NCara = Celda.ObtNCara(); double aP = 0.0; std::valarray<double> aF(NCara); TTensor<d, r + 1u> const Grad = φ.CalcGrad(Celda); TTensor<d, r> b = {}; for (std::size_t i = 0u; i < NCara; ++i) { TCara<d> const &Cara = Celda.ObtCara(i); auto const [aP_, aF_, b_] = TExprLapBase<d, r>::CalcCoef(Cara, Grad); double const Γf = Γ.CalcVal(Cara); aP += Γf * aP_; aF[i] = Γf * aF_; b += Γf * b_; } return {aP, aF, b}; } // ================================================================================================= // ============================================================================ TExprLap<d, r, void> template<std::size_t d, std::size_t r> TCoef<d, r> TExprLap<d, r, void>::CalcCoef(TCelda<d> const &Celda) const { std::size_t const NCara = Celda.ObtNCara(); double aP = 0.0; std::valarray<double> aF(NCara); TTensor<d, r + 1u> const Grad = φ.CalcGrad(Celda); TTensor<d, r> b = {}; for (std::size_t i = 0u; i < NCara; ++i) { auto const [aP_, aF_, b_] = TExprLapBase<d, r>::CalcCoef(Celda.ObtCara(i), Grad); aP += aP_; aF[i] = aF_; b += b_; } return {aP, aF, b}; } // ================================================================================================= // ========================================================================================= TExprSp template<std::size_t d, std::size_t r, typename T> TCoef<d, r> TExprSp<d, r, T>::CalcCoef(TCelda<d> const &Celda) const { double aP = 0.0; TTensor<d, r> b = {}; if (double const ρP = ρ.CalcVal(Celda); ρP >= 0.0) aP = ρP; else b = ρP * φ.CalcVal(Celda); return {aP, std::valarray<double>(Celda.ObtNCara()), b}; } // ============================================================================== FUNCIONES EN LÍNEA // =============================================================================================== d template<std::size_t d1, std::size_t r> TExprd<d1, r> inline d(TCampo<d1, r> const &φ) { return {φ}; } // ================================================================================================= // ============================================================================================= div template<std::size_t d, std::size_t r> TExprDiv<d, r - 1u> inline div(TCampo<d, r> const &φ) { return {φ}; } // ================================================================================================= template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprDiv<d, r - 1u, T<d, 1u, TRest...>> inline div(TExprBinaria<d, r, T<d, 1u, TRest...>, TCampo<d, r - 1u>, std::multiplies<>> const &Expr) { return {Expr}; } // ================================================================================================= template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprLap<d, r - 1u, T<d, 0u, TRest...>> inline div(TExprBinaria<d, r, T<d, 0u, TRest...>, TExprGrad<d, r>, std::multiplies<>> const &Expr) { return {Expr}; } // ================================================================================================= // ============================================================================================ grad template<std::size_t d, std::size_t r> TExprGrad<d, r + 1u> inline grad(TCampo<d, r> const &φ) { return {φ}; } // ================================================================================================= // =========================================================================================== gradT template<std::size_t d> TExprGradT<d, 2u> inline gradT(TCampo<d, 1u> const &φ) { return {φ}; } // ================================================================================================= // ============================================================================================= lap template<std::size_t d, std::size_t r> TExprLap<d, r> inline lap(TCampo<d, r> const &φ) { return {φ}; } // ================================================================================================= // ============================================================================================ lapT template<std::size_t d> TExprLapT<d, 1u> inline lapT(TCampo<d, 1u> const &φ) { return {φ}; } // ================================================================================================= // ============================================================================================= rot template<std::size_t d> TExprRot<d, 1u> inline rot(TCampo<d, 1u> const &φ) { return {φ}; } // ================================================================================================= // ============================================================================================== Sp template<std::size_t d, std::size_t r> TExprSp<d, r> inline Sp(TCampo<d, r> const &φ) { return {φ}; } // ================================================================================================= template<std::size_t d, std::size_t r, typename... TRest, template<std::size_t, std::size_t, typename...> typename T> TExprSp<d, r, T<d, 0u, TRest...>> Sp(TExprBinaria<d, r, T<d, 0u, TRest...>, TCampo<d, r>, std::multiplies<>> const &Expr) { return {Expr}; } } // VF namespace VF { // =========================================================================== DECLARACIÓN DE CLASES // ======================================================================================== TSistema template<std::size_t d, std::size_t r> class TSistema { private: template<std::size_t = d, std::size_t = r> class TExprΣaF : public TExprBase<TExprΣaF<d, r>> { private: TSistema const &Sistema; TCampo<d, r> const &φ; public: TExprΣaF(TSistema const &Sistema_, TCampo<d, r> const &φ_) : Sistema(Sistema_), φ(φ_) {} TTensor<d, r> CalcVal(TCelda<d> const &Celda) const { return Sistema.ΣaF(Celda.ObtID(), φ); } TCoef<d, r> CalcCoef(TCelda<d> const &) const = delete; template<std::true_type = {}> TTensor<d, r> operator [](std::size_t const i) const { return Sistema.ΣaF(i, φ); } }; // ------------------------------------------------------------------------------------------- Datos private: static inline double ε = 1e-9; // Error relativo a ||b|| static inline std::size_t MaxIt = 10'000u; // Número máximo de iteraciones std::vector<std::valarray<double>> aFVec; // Términos no diagonales public: TCampo<d, 0u> aP; // Diagonal TCampo<d, r> b; // Términos independientes // --------------------------------------------------------------------------------------- Funciones private: double static Dot(TCampo<d, r> const &lhs, TCampo<d, r> const &rhs) requires (r == 0u) { return sum(lhs * rhs); } double static Dot(TCampo<d, r> const &lhs, TCampo<d, r> const &rhs) requires (r == 1u) { return sum(lhs & rhs); } double static Dot(TCampo<d, r> const &lhs, TCampo<d, r> const &rhs) requires (r == 2u) { return sum(lhs && rhs); } TTensor<d, r> ΣaF(std::size_t const, TCampo<d, r> const &) const; void Solve(TCampo<d, r> &, double const) const; public: TSistema() = delete; template<typename T, typename U> TSistema(TExprBase<T> const &, U const &); void static DefMaxIt(std::size_t const MaxIt_) { MaxIt = MaxIt_; } std::size_t static ObtMaxIt() { return MaxIt; } void static Defε(double const ε_) { ε = ε_; } double static Obtε() { return ε; } void DefRef(std::size_t const ID, TTensor<d, r> const &Ref) { b[ID] += aP[ID] * Ref; aP[ID] *= 2.0; } TExprΣaF<d, r> ΣaF(TCampo<d, r> const &φ) const { return {*this, φ}; } // ------------------------------------------------------------------------------------------ Amigos friend class TExprΣaF<d, r>; void friend solve(TSistema const &Sistema, TCampo<d, r> &Campo, double const f = 1.0) { Sistema.Solve(Campo, f); } }; // ======================================================================== IMPLEMENTACIÓN DE CLASES // ======================================================================================== TSistema template<std::size_t d, std::size_t r> template<typename T, typename U> TSistema<d, r>::TSistema(TExprBase<T> const &lhs, U const &rhs) : aFVec(TMalla<d>::ObtNCelda()) { #pragma omp parallel for for (std::size_t i = 0u; i < TMalla<d>::ObtNCelda(); ++i) { TCelda<d> const &Celda = TMalla<d>::ObtCelda(i); TCoef const Coef = lhs.CalcCoef(Celda); aP[i] = Coef.aP; aFVec[i] = Coef.aF; if constexpr (EsExpr<U>) b[i] = rhs.CalcVal(Celda) - Coef.b; else b[i] = rhs - Coef.b; } } // ======================================================================== IMPLEMENTACIÓN DE CLASES // ======================================================================================== TSistema template<std::size_t d, std::size_t r> TTensor<d, r> TSistema<d, r>::ΣaF(std::size_t const i, TCampo<d, r> const &φ) const { TTensor<d, r> aFφ = {}; for (auto &&[ID, aF] : ranges::views::zip(TMalla<d>::ObtIDVec(i), aFVec[i])) aFφ += aF * φ[ID]; return aFφ; } // ================================================================================================= template<std::size_t d, std::size_t r> void TSistema<d, r>::Solve(TCampo<d, r> &φ, double const f) const { double const tol2 = pow<2u>(ε) * Dot(b, b); TCampo<d, r> x = φ, e = b - aP * x - ΣaF(x), e0 = e, p, s, y, z, v, t; double ρ = 1.0, α = 1.0, β, ω = 1.0, e0e0 = Dot(e0, e0); for (std::size_t It = 0u; It < MaxIt; ++It) { double const ρOld = ρ; if (Dot(e, e) < tol2) break; ρ = Dot(e0, e); if (std::abs(ρ) < pow<2u>(std::numeric_limits<double>::epsilon()) * e0e0) { e0 = e = b + aP * (1.0 - f) / f * φ - aP / f * x - ΣaF(x); ρ = e0e0 = Dot(e0, e0); } β = (ρ / ρOld) * (α / ω); p = e + β * (p - ω * v); y = p / aP; v = aP * y / f + ΣaF(y); α = ρ / Dot(e0, v); s = e - α * v; z = s / aP; t = aP * z / f + ΣaF(z); ω = Dot(t, s) / Dot(t, t); x += α * y + ω * z; e = s - ω * t; } φ = std::move(x); } } // VF namespace VF { // ============================================================================== VARIABLES EN LÍNEA // ================================================================================================= inline constexpr double C1 = 1.44; inline constexpr double C2 = 1.92; inline constexpr double σk = 1.0; inline constexpr double σε = 1.3; inline constexpr double Cμ = 0.09; inline constexpr double κ = 0.41; inline constexpr double E = 9.8; // =========================================================================== DECLARACIÓN DE CLASES // ====================================================================================== TkWallFunc template<std::size_t d, std::size_t r> requires (r == 0u) class TkWallFunc : public TNeumann<d, 0u> { public: TkWallFunc() = default; }; // ================================================================================================= // =================================================================================== TWallFuncBase template<std::size_t d> class TWallFuncBase { private: TCampo<d, 0u> const &k; // --------------------------------------------------------------------------------------- Funciones protected: TWallFuncBase(TCampo<d, 0u> const &k_) : k(k_) {} double CalcUτ(TCara<d> const &Cara) const { return std::sqrt(std::sqrt(Cμ) * k.CalcVal(Cara.ObtCeldaP())); } }; // ================================================================================================= // ====================================================================================== TεWallFunc template<std::size_t d, std::size_t r> requires (r == 0u) class TεWallFunc : public TWallFuncBase<d>, public TCCBase<d, 0u> { private: using TWallFuncBase<d>::CalcUτ; public: TεWallFunc() = delete; TεWallFunc(TCampo<d, 0u> const &k_) : TWallFuncBase<d>(k_) {} std::tuple<double, TTensor<d, 0u>> virtual CalcVal(TCara<d> const &Cara) const override { return {0.0, pow<3u>(CalcUτ(Cara)) / (κ * mag(Cara.ObtL()))}; } std::tuple<double, TTensor<d, 0u>> virtual CalcGrad(TCara<d> const &Cara) const override { return {0.0, pow<3u>(CalcUτ(Cara)) / (κ * (Cara.ObtL() & Cara.ObtL()))}; } }; // ================================================================================================= // ==================================================================================== TμEfWallFunc template<std::size_t d, std::size_t r> requires (r == 0u) class TμEfWallFunc : public TWallFuncBase<d>, public TCCBase<d, 0u> { private: double const μ; // --------------------------------------------------------------------------------------- Funciones private: using TWallFuncBase<d>::CalcUτ; public: TμEfWallFunc() = delete; TμEfWallFunc(TCampo<d, 0u> const &k_, double const μ_) : TWallFuncBase<d>(k_), μ(μ_) {} std::tuple<double, TTensor<d, 0u>> virtual CalcVal(TCara<d> const &) const override; }; // ======================================================================== IMPLEMENTACIÓN DE CLASES // ==================================================================================== TμEfWallFunc template<std::size_t d, std::size_t r> std::tuple<double, TTensor<d, 0u>> TμEfWallFunc<d, r>::CalcVal(TCara<d> const &Cara) const { double const yPlus = CalcUτ(Cara) * mag(Cara.ObtL()) / μ; return {0.0, μ * κ * yPlus / std::log(E * yPlus)}; } } // VF using namespace VF; using namespace std; // -------------------------------------------------------------------------------------- Constantes constexpr TVector<2u> U0 = {10.0, 0.0}; constexpr double k0 = 0.375, ε0 = 14.855; constexpr double ν = 1.8e-5; constexpr double α = 0.7; int main() { // ------------------------------------------------------------------------------------------- Malla TMalla2D::ReadMSH("backwardStep.msh"); // ------------------------------------------------------------------------------------------ Campos TCampoVectorial2D U; TCampoEscalar2D p, k, ε, νt, νEf; // ------------------------------------------------------------------------- Condiciones de contorno U.DefCC<TDirichlet>("inlet", U0); k.DefCC<TDirichlet>("inlet", k0); ε.DefCC<TDirichlet>("inlet", ε0); U.DefCC<TDirichlet>("wall"); k.DefCC<TkWallFunc>("wall"); ε.DefCC<TεWallFunc>("wall", k); νEf.DefCC<TμEfWallFunc>("wall", k, ν); p.DefCC<TDirichlet>("outlet"); // --------------------------------------------------------------------------- Condiciones iniciales U = U0; k = k0; ε = ε0; while (true) { // ---------------------------------------------------------------------------------- Momento lineal νt = Cμ * k * k / ε; νEf = ν + νt; TSistema const UEc = div(U * U) - div(νEf * grad(U)) == (gradT(U) & grad(νt)) - grad(p); solve(UEc, U, α); // ------------------------------------------------------------------------------------- Continuidad U = TCampo((UEc.b + grad(p) - UEc.ΣaF(U)) / UEc.aP); TCampo const pOld = p; solve(div((1.0 / UEc.aP) * grad(p)) == div(U), p); U -= grad(p) / UEc.aP; if (sum(mag(p - pOld)) < 0.1) break; p = lerp(p, pOld, α); // ----------------------------------------------------------------------------- Transporte de k y ε TCampo const ω = ε / k; auto const G = νt * ((grad(U) + gradT(U)) && grad(U)); // No se evalúa la expresión solve(div(U * k) - div(νt * grad(k)) / σk + ω * Sp(k) == G, k, α); solve(div(U * ε) - div(νt * grad(ε)) / σε + C2 * ω * Sp(ε) == C1 * ω * G, ε, α); } return 0; }
Become a Patron
Sponsor on GitHub
Donate via PayPal
Compiler Explorer Shop
Source on GitHub
Mailing list
Installed libraries
Wiki
Report an issue
How it works
Contact the author
CE on Mastodon
CE on Bluesky
Statistics
Changelog
Version tree