FEENOX(1)

๐Ÿ“ฆ
version 1.2 ยท 2 versions
, in ubuntu resolute, ubuntu stonking

๐Ÿ’– FEENOX ๐Ÿ’–

okay bestie, sit down, hydrate, because this is the FeenoX User Manual and it's literally the whole vibe. ๐Ÿ’…


โœจ NAME โœจ

FeenoX โ€” a cloud-first free no-X uniX-like finite-element(ish) computational engineering tool ๐Ÿฆ‹

Like. it's giving "I solve your engineering problems and I don't even need a GUI". periodt.


๐ŸŒธ SYNOPSIS ๐ŸŒธ

the basic usage is to execute the feenox binary passing a path to an input file that defines the problem, along with other options and command-line replacement arguments which are explained below:

feenox [options ...]  input-file
[optional_commandline_replacement_arguments ...]

and for large problems that do not fit in a single computer, a parallel run using mpirun(1) will be needed:

mpirun -n number_of_threads feenox [options ...]  input-file
[optional_commandline_replacement_arguments ...]

๐Ÿ’œ DESCRIPTION ๐Ÿ’œ

so. FeenoX is a computational tool that can solve engineering problems which are usually casted as differential-algebraic equations (DAEs) or partial differential equations (PDEs). and the way I like to say it is: it's to finite elements programs and libraries what Markdown is to Word and TeX, respectively. iconic behavior fr fr. in particular, it can solve:

  • dynamical systems defined by a set of user-provided DAEs (such as plant control dynamics for example) ๐ŸŽ€
  • mechanical elasticity
  • heat conduction
  • structural modal analysis
  • neutron diffusion
  • neutron transport

FeenoX reads a plain-text input file which contains the problem definition and writes 100%-user defined results in ASCII (through PRINT or other user-defined output instructions within the input file). for PDE problems, it needs a reference to at least one Gmsh mesh file for the discretization of the domain. and it can write post-processing views in either .msh, .vtu or .vtk formats. no cap.

keep in mind that FeenoX is just a back end reading a set of input files and writing a set of output files following the design philosophy of Unix (separation, composition, representation, economy, extensibility, etc). think of it as a transfer function (or a filter in computer-science jargon) between input files and output files:

                                    +------------+
        mesh (*.msh)  }             |            |             { terminal
        data (*.dat)  } input ----> |   FeenoX   |----> output { data files
        input (*.fee) }             |            |             { post (vtk/msh)
                                    +------------+

following the Unix programming philosophy, there are no graphical interfaces attached to the FeenoX core, although a wide variety of pre and post-processors can be used with FeenoX. to illustrate the transfer-function approach, consider the following input file that solves Laplace's equation โˆ‡^2^ฯ• = 0 on a square with some space-dependent boundary conditions:

ฯ•(x, y) = +y   for x = โˆ’1 (left)

ฯ•(x, y) = โˆ’y   for x = +1 (right)

โˆ‡ฯ• โ‹… nฬ‚ = sin (ฯ€/2x)   for y = โˆ’1 (bottom)

โˆ‡ฯ• โ‹… nฬ‚ = 0   for y = +1 (top)
       PROBLEM laplace 2d
       READ_MESH square-centered.msh # [-1:+1]x[-1:+1]

       # boundary conditions
       BC left    phi=+y
       BC right   phi=-y
       BC bottom  dphidn=sin(pi/2*x)
       BC top     dphidn=0

       SOLVE_PROBLEM

       # same output in .msh and in .vtk formats
       WRITE_MESH laplace-square.msh phi VECTOR dphidx dphidy 0
       WRITE_MESH laplace-square.vtk phi VECTOR dphidx dphidy 0

Laplace's equation solved with FeenoX ๐Ÿ˜Œ

the .msh file can be post-processed with Gmsh, and the .vtu/.vtk file can be post-processed with Paraview. and see https://www.caeplex.com for a mobile-friendly web-based interface for solving finite elements in the cloud directly from the browser. cloud-first, we been knew. โ˜๏ธโœจ


โœจ OPTIONS โœจ

feenox [options] inputfile [replacement arguments] [petsc options]
  • -h, --help โ€” display options and detailed explanations of command-line usage ๐Ÿ™‹
  • -v, --version โ€” display brief version information and exit
  • -V, --versions โ€” display detailed version information
  • -c, --check โ€” validates if the input file is sane or not (vibe check for your input file, basically ๐Ÿซถ)
  • --pdes โ€” list the types of PROBLEMs that FeenoX can solve, one per line
  • --elements_info โ€” output a document with information about the supported element types
  • --ast โ€” dump an abstract syntax tree of the input
  • --linear โ€” force FeenoX to solve the PDE problem as linear
  • --non-linear โ€” force FeenoX to solve the PDE problem as non-linear
  • --progress โ€” print ASCII progress bars when solving PDEs (the little bars, the drama, obsessed)
  • --mumps โ€” ask PETSc to use the direct linear solver MUMPS
  • --gamg โ€” ask PETSc to use a GAMG-preconditioned iterative linear solver

instructions will be read from standard input if - is passed as inputfile, i.e.

$ echo 'PRINT 2+2' | feenox -
4

the optional [replacement arguments] part of the command line mean that each argument after the input file that does not start with an hyphen will be expanded verbatim in the input file in each occurrence of $1, $2, etc. for example

$ echo 'PRINT $1+$2' | feenox - 3 4
7

PETSc and SLEPc options can be passed in [petsc options] (or [options]) as well, with the difference that two hyphens have to be used instead of only once. for example, to pass the PETSc option -ksp_view the actual FeenoX invocation should be

$ feenox input.fee --ksp_view

for PETSc options that take values, en equal sign has to be used:

$ feenox input.fee --mg_levels_pc_type=sor

see https://www.seamplex.com/feenox/examples for annotated examples. ๐Ÿงš


๐ŸŒˆ EXAMPLES ๐ŸŒˆ

(nothing here in this part, bestie โ€” scroll on)


๐Ÿ’€ EXIT STATUS ๐Ÿ’€

(empty for now)


๐Ÿ“ ENVIRONMENT ๐Ÿ“

(empty for now)


๐ŸŽ€ FILES ๐ŸŽ€

(empty for now)


โญ CONFORMING TO โญ

(empty for now)


๐Ÿ‘‘ INPUT-FILE KEYWORDS ๐Ÿ‘‘

๐Ÿ’… GENERIC KEYWORDS ๐Ÿ’…

๐Ÿ›‘ ABORT

catastrophically abort the execution and quit FeenoX. yes, catastrophically. that's the word they used. ๐Ÿ˜ญ

ABORT

whenever the instruction ABORT is executed, FeenoX quits with a non-zero error leve. it does not close files nor unlock shared memory objects. the objective of this instruction is to either debug complex input files by using only parts of them or to conditionally abort the execution using IF clauses.

bestie this is the "I'm out" button. it's dramatic on purpose. ๐Ÿ”ฅ

โœจ ALIAS

define a scalar alias of an already-defined identifier.

ALIAS { <new_var_name> IS <existing_object> | <existing_object> AS <new_name> }

the existing object can be a variable, a vector element or a matrix element. in the first case, the name of the variable should be given as the existing object. in the second case, to alias the second element of vector v to the new name new, v(2) should be given as the existing object. in the third case, to alias second element (2,3) of matrix M to the new name new, M(2,3) should be given as the existing object.


๐Ÿ”’ CLOSE

explicitly close a file after input/output.

CLOSE <name>

the given <name> can be either a fixed-string path or an already-defined FILE.


๐Ÿซถ DEFAULT_ARGUMENT_VALUE

give a default value for an optional commandline argument.

DEFAULT_ARGUMENT_VALUE <constant> <string>

if a $n construction is found in the input file but the commandline argument was not given, the default behavior is to fail complaining that an extra argument has to be given in the commandline. with this keyword, a default value can be assigned if no argument is given, thus avoiding the failure and making the argument optional. the <constant> should be 1, 2, 3, etc. and <string> will be expanded character-by-character where the $n construction is. whether the resulting expression is to be interpreted as a string or as a numerical expression will depend on the context.


๐Ÿ“ FILE

define a file with a particularly formatted name to be used either as input or as output.

< FILE | OUTPUT_FILE | INPUT_FILE > <name> PATH <format> expr_1 expr_2 ... expr_n [ INPUT | OUTPUT | APPEND | MODE <fopen_mode> ]

for reading or writing into files with a fixed path, this instruction is usually not needed as the FILE keyword of other instructions (such as PRINT or MESH) can take a fixed-string path as an argument. however, if the file name changes as the execution progresses (say because one file for each step is needed), then an explicit FILE needs to be defined with this keyword and later referenced by the given name.

the path should be given as a printf-like format string followed by the expressions which should be evaluated in order to obtain the actual file path. the expressions will always be floating-point expressions, but the particular integer specifier %d is allowed and internally transformed to %.0f. the file can be explicitly defined and INPUT, OUTPUT or a certain fopen() mode can be given (i.e. โ€œaโ€). if not explicitly given, the nature of the file will be taken from context, i.e. FILEs in PRINT will be OUTPUT and FILEs in FUNCTION will be INPUT.

this keyword just defines the FILE, it does not open it. the file will be actually opened (and eventually closed) automatically. in the rare case where the automated opening and closing does not fit the expected workflow, the file can be explicitly opened or closed with the instructions FILE_OPEN and FILE_CLOSE.


๐Ÿ“ˆ FIT

find parameters to fit an analytical function to a pointwise-defined function.

FIT <function_to_be_fitted>  TO <function_with_data> VIA <var_1> <var_2> ... <var_n>
 [ GRADIENT <expr_1> <expr_2> ... <expr_n> ]
 [ RANGE_MIN <expr_1> <expr_2> ... <expr_j> ]
 [ RANGE_MAX <expr_1> <expr_2> ... <expr_n> ]
 [ TOL_REL <expr> ] [ TOL_ABS <expr> ] [ MAX_ITER <expr> ]
 [ VERBOSE ]

more parts incoming, don't go anywhere ๐Ÿ’–

okay so the function that's holding the data? It has to be point-wise defined ๐Ÿ’… โ€” like a FUNCTION that got read from a file, or one with inline DATA, or one defined over a mesh. And the function you're fitting has to be parametrized with at least one of the variables you dropped after the USING keyword, periodt.

Say you wanna fit f(x, y) = ax^2^ + bsqrt(y) to a pointwise-defined g(x, y). You'd write FIT f TO g VIA a b. Only the names of the functions, bestie โ€” not the arguments. And both functions need the same number of arguments, no exceptions. ๐Ÿซถ

The initial guess of the solution is literally just the initial value of the variables after the VIA keyword. Want to hand over analytical expressions for the gradient of the function you're fitting, with respect to the parameters being fitted? You can, with the GRADIENT keyword. โœจ If you don't? No stress, it just computes the gradient numerically using finite differences.

You can also give a range over which the residuals get minimized with RANGE_MIN and RANGE_MAX. Heads up: those expressions give the range of the arguments of the functions, not of the parameters โ€” super important. For multidimensional fits the range is a hypercube. And if you give no range at all, then all the definition points of the data function get used for the fit. ๐ŸŒธ

Convergence control? Got you: relative and absolute tolerances via TOL_REL (default DEFAULT_NLIN_FIT_EPSREL) and TOL_ABS (default DEFAULT_NLIN_FIT_EPSABS), plus the maximum number of iterations with MAX_ITER (default DEFAULT_NLIN_FIT_MAX_ITER). And if you throw in the optional VERBOSE keyword, some data from the intermediate steps gets written to standard output. ๐Ÿ‘‘


โœจ FUNCTION โœจ

Define a scalar function of one or more variables. Iconic behavior, honestly.

FUNCTION <function_name>(<var_1>[,var2,...,var_n]) {
   = <expr> |
   FILE { <file> } |
   VECTORS <vector_1> <vector_2> ... <vector_n> <vector_data> |
   MESH <mesh> |
   DATA <num_1> <num_2> ... <num_N>
 }
 [ COLUMNS <expr_1> <expr_2> ... <expr_n> <expr_n+1> ]
 [ INTERPOLATION { linear | polynomial | spline | spline_periodic | akima | akima_periodic | steffen |
 nearest | shepard | shepard_kd | bilinear } ]
 [ INTERPOLATION_THRESHOLD <expr> ] [ SHEPARD_RADIUS <expr> ] [ SHEPARD_EXPONENT <expr> ]

The number of variables n comes from the number of arguments you put between the parentheses after the function name. Those arguments get defined as new variables if they hadn't already been explicitly defined as scalar variables. Not me explaining scoping in a manpage, but here we are. ๐Ÿ™ˆ

If your function is an algebraic expression, you can use the short-hand operator = (or := if you're feeling Maxima-compatible). So like, FUNCTION f(x) = x^2 is equivalent to f(x) = x^2 (or f(x) := x^2). โœจ

  • Give it a FILE and it expects an ASCII file with at least n + 1 columns. By default the first n columns are the values of the arguments and the last column is the value of the function at those points. You can reorder that with the COLUMNS keyword, which expects n + 1 expressions corresponding to the column numbers.
  • Give it VECTORS and it expects a set of n + 1 vectors of the same size. The first n are the arguments, the last one is the function values.
  • Give it MESH and the function is point-wise defined over the mesh topology. Translation: the independent variables (i.e. the spatial coordinates) coincide with the mesh nodes, and the dependent variable (i.e. the function value) is set by โ€œfillingโ€ a vector named vec_f (where f has to be replaced with the function name) of size equal to the number of nodes. ๐Ÿฆ‹

You can also define the function pointwise inline in the input using DATA. This one should be the last keyword of the line, followed by N = k โ‹… (n + 1) expressions giving k definition points: n arguments and the value of the function. Multiline continuation using brackets { and } can be used for a clean data organization. ๐Ÿ’…

Interpolation schemes can be given for either one- or multi-dimensional functions with INTERPOLATION.

Available schemes for n = 1 are:

  • linear
  • polynomial, the grade is equal to the number of data minus one
  • spline, cubic (needs at least 3 points)
  • spline_periodic
  • akima (needs at least 5 points)
  • akima_periodic (needs at least 5 points)
  • steffen, always-monotonic splines-like interpolator

Default interpolation scheme for one-dimensional functions is DEFAULT_INTERPOLATION.

Available schemes for n > 1 are:

  • nearest, f(x) is equal to the value of the closest definition point
  • shepard, inverse distance weighted average definition points (might lead to inefficient evaluation)
  • shepard_kd, average of definition points within a kd-tree (more efficient evaluation provided SHEPARD_RADIUS is set to a proper value)
  • bilinear, only available if the definition points configure an structured hypercube-like grid. If n > 3, SIZES should be given.

For n > 1, if the euclidean distance between the arguments and the definition points is smaller than INTERPOLATION_THRESHOLD, the definition point is returned and no interpolation is performed. Default value is square root of DEFAULT_MULTIDIM_INTERPOLATION_THRESHOLD.

The initial radius of points to take into account in shepard_kd is given by SHEPARD_RADIUS. If no points are found, the radius is double until at least one definition point is found. The radius is doubled until at least one point is found. Default is DEFAULT_SHEPARD_RADIUS. The exponent of the shepard method is given by SHEPARD_EXPONENT. Default is DEFAULT_SHEPARD_EXPONENT. ๐ŸŒˆ


โœจ IF โœจ

Execute a set of instructions if a condition is met. Simple as that, bestie.

IF expr
  <block_of_instructions_if_expr_is_true>
 [ ELSE
  <block_of_instructions_if_expr_is_false> ]
 ENDIF

โœจ IMPLICIT โœจ

Define whether implicit definition of variables is allowed or not.

IMPLICIT { NONE | ALLOWED }

By default, FeenoX allows variables (but not vectors nor matrices) to be implicitly declared. But likeโ€ฆ typos happen, and typos become bugs, so you can force explicit declaration of variables by giving IMPLICIT NONE. Whether implicit declaration is allowed or explicit declaration is required depends on the last IMPLICIT keyword given, which by default is ALLOWED. ๐ŸŽ€


โœจ INCLUDE โœจ

Include another FeenoX input file. It's giving copy-paste but make it fancy.

INCLUDE <file_path> [ FROM <num_expr> ] [ TO <num_expr> ]

Includes the input file located in the string file_path at the current location. The effect is the same as copying and pasting the contents of the included file at the location of the INCLUDE keyword. The path can be relative or absolute. ๐Ÿ’ซ

Note, however, that when including files inside IF blocks that instructions are conditionally-executed but all definitions (such as function definitions) are processed at parse-time independently from the evaluation of the conditional. Read that twice, bestie. ๐Ÿ˜ญ

The included file has to be an actual file path (i.e. it cannot be a FeenoX FILE) because it needs to be resolved at parse time. Yet, the name can contain a commandline replacement argument such as $1 so INCLUDE $1.fee will include the file specified after the main input file in the command line. The optional FROM and TO keywords can be used to include only portions of a file. ๐Ÿ‘‘


โœจ MATRIX โœจ

Define a matrix.

MATRIX <name> ROWS <expr> COLS <expr> [ DATA <expr_1> <expr_2> ... <expr_n> |

A new matrix of the prescribed size is defined. The number of rows and columns can be an expression which will be evaluated the very first time the matrix is used and then kept at those constant values. All elements will be initialized to zero unless DATA is given (which should be the last keyword of the line), in which case the expressions will be evaluated the very first time the matrix is used and row-major-assigned to each of the elements. If there are less elements than the matrix size, the remaining values will be zero. If there are more elements than the matrix size, the values will be ignored. No cap. ๐Ÿ”ฅ


โœจ OPEN โœจ

Explicitly open a file for input/output.

OPEN <name> [ MODE <fopen_mode> ]

The given <name> can be either a fixed-string path or an already-defined FILE. The mode is only taken into account if the file is not already defined. Default is write w. ๐Ÿ’…


โœจ PRINT โœจ

Write plain-text and/or formatted data to the standard output or into an output file. Main character energy for your data, honestly.

PRINT [ <object_1> <object_2> ... <object_n> ] [ TEXT <string_1> ... TEXT <string_n> ]
 [ FILE { <file_path> | <file_id> } ] [ HEADER ] [ NONEWLINE ] [ SEP <string> ]
 [ SKIP_STEP <expr> ] [ SKIP_STATIC_STEP <expr> ] [ SKIP_TIME <expr> ] [ SKIP_HEADER_STEP <expr> ]

okay so every argument object that is not a keyword of the PRINT instruction? bestie, it's going into the output. no notes. ๐Ÿ’…

Those objects can be a matrix, a vector, or literally any valid scalar algebraic expression. And if the thing you handed over can't be solved into a valid matrix, vector or expression, it gets treated as a string literal if IMPLICIT is ALLOWED โ€” otherwise the parser just raises an error at you. ๐Ÿซ 

Want to force a literal string even when it does resolve to a totally valid numerical expression? Slap the TEXT keyword in front of it. Like PRINT TEXT 1+1 prints 1+1 instead of 2. iconic, honestly. Objects and string literals can be mixed together and given in any order. no rules, we're free. ๐Ÿฆ‹

Now โ€” hashes # showing up literally inside text strings have to be quoted, otherwise the parser treats them as comments within the FeenoX input file and just ignores the rest of the line. so it's PRINT "\# this is a printed comment". iykyk. ๐ŸŽ€

Whenever an argument starts with a percentage sign %, it's treated as a C printf-compatible format specifier, and all the objects that follow it get printed using that format until a new format definition is found. The objects are treated as double-precision floating point numbers, so bestie, only give floating point formats. See the printf(3) man page for further details. The default format is DEFAULT_PRINT_FORMAT.

Matrices, vectors, scalar expressions, format modifiers and string literals can be given in any desired order, and they're processed from left to right. Vectors are printed element-by-element in a single row. (Peep PRINT_VECTOR if you want one or more vectors printed with one element per line, i.e. vertically.) Matrices are printed element-by-element in a single line using row-major ordering if they're mixed with other objects โ€” but in the natural row and column fashion if the matrix is the only given object in the PRINT instruction. it's giving โœจrangeโœจ

If the FILE keyword is not provided, the default is to write to stdout. If the HEADER keyword is given, a single line containing the literal text given for each object is printed at the very first time the PRINT instruction is processed, starting with a hash # character.

If the NONEWLINE keyword is not provided, the default is to write a newline \n character after all the objects are processed. Otherwise, if the last token to be printed is a numerical value, a separator string will be printed but not the newline \n character. If the last token is a string, neither the separator nor the newline will be printed. The SEP keyword expects a string used to separate printed objects. To print objects without any separation in between, give an empty string like SEP "". The default is a tabulator character DEFAULT_PRINT_SEPARATOR character. To print an empty line, write PRINT with no arguments at all. ๐Ÿซถ

By default the PRINT instruction is evaluated every step. If the SKIP_STEP (SKIP_STATIC_STEP) keyword is given, the instruction is processed only every the number of transient (static) steps that results in evaluating the expression, which may not be constant. The SKIP_HEADER_STEP keyword works similarly for the optional HEADER, but by default it is only printed once. The SKIP_TIME keyword uses time advancements to choose how to skip printing and may be useful for non-constant time-step problems. ๐Ÿ‘‘


โœจ PRINTF โœจ

Instruction akin to C's printf. Instruction akin to C's printf executed locally from all MPI ranks. so like, everyone's doing it, but locally. ๐Ÿ˜Œ

PRINTF PRINTF_ALL format_string [ expr_1 [ expr_2 [ ... ] ] ]

The format_string should be a printf-like string containing double-precision format specifiers. A matching number of expressions should be given. No newline is written if not explicitly asked for in the format string with \n.

Do not ask for string literals %s. ๐Ÿšซ

As always, to get a literal % use %% in the format string.


Print one or more functions as a table of values of dependent and independent variables. it's giving "data table, but make it cute." ๐Ÿ“Š

PRINT_FUNCTION <function_1> [ { function | expr } ... { function | expr } ]
 [ FILE { <file_path> | <file_id> } ] [ HEADER ]
 [ MIN <expr_1> <expr_2> ... <expr_k> ] [ MAX <expr_1> <expr_2> ... <expr_k> ]
 [ STEP <expr_1> <expr_2> ... <expr_k> ] [ NSTEPs <expr_1> <expr_2> ... <expr_k> ]
 [ FORMAT <print_format> ] <vector_1> [ { vector | expr } ... { vector | expr } ]

Each argument should be either a function or an expression. The output of this instruction consists of n + k columns, where n is the number of arguments of the first function of the list and k is the number of functions and expressions given. The first n columns are the arguments (independent variables) and the last k one has the evaluated functions and expressions. The columns are separated by a tabulator, which is the format that most plotting tools understand. no cap. ๐Ÿงš

Only function names without arguments are expected. All functions should have the same number of arguments. Expressions can involve the arguments of the first function.

If the FILE keyword is not provided, the default is to write to stdout. If HEADER is given, the output is prepended with a single line containing the names of the arguments and the names of the functions, separated by tabs. The header starts with a hash # that usually acts as a comment and is ignored by most plotting tools.

If there is no explicit range where to evaluate the functions and the first function is point-wise defined, they are evaluated at the points of definition of the first one. The range can be explicitly given as a product of n ranges [x_i,min, x_i,max] for i = 1, ..., n. The values x_i,min and x_i,max are given with the MIN and MAX keywords. The discretization steps of the ranges are given by either STEP that gives ฮดx or NSTEPS that gives the number of steps. If the first function is not point-wise defined, the ranges are mandatory. ๐Ÿซก


Print the elements of one or more vectors, one element per line. vertical vibes only. โฌ‡๏ธ

PRINT_VECTOR
 [ FILE { <file_path> | <file_id> } ] [ HEADER ]
 [ SEP <string> ]

Each argument should be either a vector or an expression of the integer i. If the FILE keyword is not provided, the default is to write to stdout. If HEADER is given, the output is prepended with a single line containing the names of the arguments and the names of the functions, separated by tabs. The header starts with a hash # that usually acts as a comment and is ignored by most plotting tools.

The SEP keyword expects a string used to separate printed objects. To print objects without any separation in between, give an empty string like SEP "". The default is a tabulator character DEFAULT_PRINT_SEPARATOR character.


๐Ÿ’– READ_DATA ๐Ÿ’–

Read data from a file (or standard output, or pipe, or queue, etc.). she's flexible, she gets around. ๐Ÿ’โ€โ™€๏ธ

READ_DATA FILE { <file_path> | <file_id> } [ <variable1> [ <variable2> [...] [ <vector1> [...] ] ] ]

๐Ÿ”ฅ SOLVE ๐Ÿ”ฅ

Solve a (small) system of non-linear equations. small is doing a lot of work in that sentence, bestie.

SOLVE FOR <n> UNKNOWNS <var_1> <var_2> ... <var_n> [ METHOD { dnewton | hybrid | hybrids | broyden } ]
 [ EPSABS <expr> ] [ EPSREL <expr> ] [ MAX_ITER <expr> ]

๐ŸŽ€ SORT_VECTOR ๐ŸŽ€

Sort the elements of a vector, optionally making the same rearrangement in another vector. chaos โ†’ order. โœจ

SORT_VECTOR <vector> [ ASCENDING | DESCENDING ] [ <other_vector> ]

This instruction sorts the elements of <vector> into either ascending or descending numerical order. If <other_vector> is given, the same rearrangement is made on it. Default is ascending order.


๐Ÿ’œ VAR ๐Ÿ’œ

Explicitly define one or more scalar variables. sometimes you just gotta declare yourself. ๐Ÿ’…

VAR <name_1> [ <name_2> ] ... [ <name_n> ]

When implicit definition is allowed (see IMPLICIT), scalar variables need not to be defined before being used if from the context FeenoX can tell that a scalar variable is needed. For instance, when defining a function like f(x) = x^2 it is not needed to declare x explicitly as a scalar variable. But if one wants to define a function like g(x) = integral(f(x'), x', 0, x) then the variable x' needs to be explicitly defined as VAR x' before the integral.


๐ŸŒˆ VECTOR ๐ŸŒˆ

Define a vector. main character energy, in vector form. ๐Ÿ‘‘

VECTOR <name> SIZE <expr> [ FUNCTION_DATA <function> ] [ DATA <expr_1> <expr_2> ... <expr_n> |

A new vector of the prescribed size is defined. The size can be an expression which will be evaluated the very first time the vector is used and then kept at that constant value. If the keyword FUNCTION_DATA is given, the elements of the vector will be synchronized with the independent values of the function, which should be point-wise defined. The sizes of both the function and the vector should match.

All elements will be initialized to zero unless DATA is given (which should be the last keyword of the line), in which case the expressions will be evaluated the very first time the vector is used and assigned to each of the elements. If there are less elements than the vector size, the remaining values will be zero. If there are more elements than the vector size, the values will be ignored. ๐Ÿ’€


๐Ÿ“ INITIAL_CONDITIONS ๐Ÿ“

Define how initial conditions of DAE problems are computed.

INITIAL_CONDITIONS { AS_PROVIDED | FROM_VARIABLES | FROM_DERIVATIVES }

In DAE problems, initial conditions may be either:

  • equal to the provided expressions (AS_PROVIDED)
  • the derivatives computed from the provided phase-space variables (FROM_VARIABLES)
  • the phase-space variables computed from the provided derivatives (FROM_DERIVATIVES)

In the first case, it is up to the user to fulfill the DAE system at t = 0. If the residuals are not small enough, a convergence error will occur. (so like, don't do that, or do it correctly. ๐Ÿ™ˆ) The FROM_VARIABLES option means calling IDA's IDACalcIC routine with the parameter IDA_YA_YDP_INIT. The FROM_DERIVATIVES option means calling IDA's IDACalcIC routine with the parameter IDA_Y_INIT.

Wasora should be able to automatically detect which variables in phase-space are differential and which are purely algebraic. However, the DIFFERENTIAL keyword may be used to explicitly define them. See the SUNDIALS documentation for further information.

โญ PHASE_SPACE โญ

Ask FeenoX to solve a set of algebraic-differential equations and define the variables, vectors and/or matrices that span the phase space.

PHASE_SPACE PHASE_SPACE <vars> ... <vectors> ... <matrices> ...

๐Ÿงš TIME_PATH ๐Ÿงš

Force time-dependent problems to pass through specific instants of time. like a lil itinerary for your simulation. ๐Ÿ—บ๏ธ

TIME_PATH <expr_1> [ <expr_2>  [ ... <expr_n> ] ]

The time step dt will be reduced whenever the distance between the current time t and the next expression in the list is greater than dt so as to force t to coincide with the expressions given. The list of expressions should evaluate to a sorted list of values for all times.


๐ŸŒธ BC ๐ŸŒธ

Define a boundary condition to be applied to faces, edges and/or vertices. setting boundaries, literally. ๐Ÿซถ

BC <name> [ MESH <name> ] [ GROUP <name_1>  GROUP <name_2> ... ] [ <bc_data1> <bc_data2> ... ] [ GROUPS <name_1>  <name_2> ... ]

If the name of the boundary condition matches a physical group in the mesh, and you didn't pass GROUP or GROUPS, bestie it just links itself up automatically to that physical group. No manual labor required. ๐Ÿ’…

If there are many meshes floating around, the mesh this keyword refers to has to be given with MESH.

And if the boundary condition applies to more than one physical group in the mesh? Add as many GROUP keywords as you need. Stack 'em. ๐Ÿ‘‘

Each <bc_data> argument is a single string whose meaning depends on the type of problem being solved. For instance T=150*sin(x/pi) prescribes the temperature to depend on space as the provided expression in a thermal problem, and fixed fixes the displacements in all the directions in a mechanical or modal problem. See the particular section on boundary conditions for further details. ๐Ÿฆ‹

If the keyword GROUPS is given, then the rest of the tokens are parsed as group names where the boundary condition is applied. And heads up โ€” if either GROUP or GROUPS are given explicitly, then the BC name is not used to try to implicitly link it to a physical group in the mesh. The name stays in its lane. ๐Ÿซถ


๐Ÿ’ฅ COMPUTE_REACTION ๐Ÿ’ฅ

Compute the reaction (force, moment, power, etc.) at selected face, edge or vertex. It's giving "how hard is this thing pushing back" energy. ๐Ÿ”ฅ

COMPUTE_REACTION <physical_group> [ MOMENT [ X0 <expr> ] [ Y0 <expr> ] [ Z0 <expr> ] ] RESULT { <variable> | <vector> }

If the MOMENT keyword is not given, the zero-th order reaction is computed โ€” i.e. force in elasticity and power in thermal. If MOMENT is given, then the coordinates of the center can be given with X0, Y0 and Z0. And if you don't hand those over? The moment is computed about the barycenter of the physical group. It'll figure it out itself, main character behavior. โœจ

The resulting reaction will be stored in the variable (thermal) or vector (elasticity) provided. If the variable or vector does not exist, it will be created. Yes queen, manifest it. ๐ŸŒธ


๐Ÿ—‚๏ธ DUMP ๐Ÿ—‚๏ธ

Dump raw PETSc objects used to solve PDEs into files. Raw. Unfiltered. No cap. ๐Ÿ’€

DUMP [ FORMAT { binary | ascii | octave } ] [ K |   K_bc |   b |   b_bc |   M |   M_bc |

๐Ÿ” FIND_EXTREMA ๐Ÿ”

Find and/or compute the absolute extrema of a function or expression over a mesh (or a subset of it). Peak detection but make it math. โ›ฐ๏ธ

FIND_EXTREMA { <expression> | <function> } [ OVER <physical_group> ] [ MESH <mesh_identifier> ] [ NODES | CELLS | GAUSS ]
  [ MIN <variable> ] [ MAX <variable> ] [ X_MIN <variable> ] [ X_MAX <variable> ] [ Y_MIN <variable> ] [ Y_MAX <variable> ] [ Z_MIN <variable> ] [ Z_MAX <variable> ] [ I_MIN <variable> ] [ I_MAX <variable> ]

Either an expression or a function of space x, y and/or z should be given. By default the search is performed over the highest-dimensional elements of the mesh โ€” i.e. over the whole volume, area or length for three, two and one-dimensional meshes, respectively.

If the search is to be carried out over just a physical group, it has to be given in OVER. If there are more than one mesh defined, an explicit one has to be given with MESH.

If neither NODES, CELLS or GAUSS is given then the search is performed over the three of them. โœจ

  • With NODES only the function or expression is evaluated at the mesh nodes.
  • With CELLS only the function or expression is evaluated at the element centers.
  • With GAUSS only the function or expression is evaluated at the Gauss points.

The value of the absolute minimum (maximum) is stored in the variable indicated by MIN (MAX). If the variable does not exist, it is created.

The value of the x-y-z coordinate of the absolute minimum (maximum) is stored in the variable indicated by X_MIN-Y_MIN-Z_MIN (X_MAX-Y_MAX-Z_MAX). If the variable does not exist, it is created.

The index (either node or cell) where the absolute minimum (maximum) is found is stored in the variable indicated by I_MIN (I_MAX). ๐Ÿ“


๐Ÿงฎ INTEGRATE ๐Ÿงฎ

Spatially integrate a function or expression over a mesh (or a subset of it). Math girlies rise. ๐Ÿ’œ

INTEGRATE { <expression> | <function> } [ OVER <physical_group> ] [ MESH <mesh_identifier> ] [ GAUSS | CELLS ]
  RESULT <variable>

Either an expression or a function of space x, y and/or z should be given. If the integrand is a function, do not include the arguments โ€” i.e. instead of f(x,y,z) just write f. The results should be the same but efficiency will be different (faster for pure functions).

By default the integration is performed over the highest-dimensional elements of the mesh, i.e. over the whole volume, area or length for three, two and one-dimensional meshes, respectively. If the integration is to be carried out over just a physical group, it has to be given in OVER. If there are more than one mesh defined, an explicit one has to be given with MESH.

Either GAUSS or CELLS define how the integration is to be performed:

  • With GAUSS the integration is performed using the Gauss points and weights associated to each element type.
  • With CELLS the integral is computed as the sum of the product of the integrand at the center of each cell (element) and the cell's volume.
Do expect differences in the results and efficiency between these two approaches depending on the nature of the integrand. It's not a bug, it's a lifestyle. ๐Ÿ™ˆ

The scalar result of the integration is stored in the variable given by the mandatory keyword RESULT. If the variable does not exist, it is created. ๐ŸŽ€


๐Ÿงฑ MATERIAL ๐Ÿงฑ

Define a material and its properties to be used in volumes.

MATERIAL <name> [ MESH <name> ] [ LABEL <name_1>  [ LABEL <name_2> [ ... ] ] ]
 [ MODEL <model_name> ]
 [ <property_name_1>=<expr_1> [ <property_name_2>=<expr_2> [ ... ] ] ]

If the name of the material matches a physical group in the mesh, it is automatically linked to that physical group. If there are many meshes, the mesh this keyword refers to has to be given with MESH. If the material applies to more than one physical group in the mesh, they can be added using as many LABEL keywords as needed. โœจ

For the mechanical problem, a string with the material model can be passed with MODEL. See the particular documentation of the mechanical problem for details.

The names of the properties in principle can be arbitrary, but each problem type needs a minimum set of properties defined with particular names. For example, steady-state thermal problems need at least the conductivity which should be named k. If the problem is transient, it will also need heat capacity rhocp or diffusivity alpha. Mechanical problems need Young modulus E and Poisson's ratio nu. Modal also needs density rho. Check the particular documentation for each problem type โ€” iykyk. ๐Ÿ“š

Besides these mandatory properties, any other one can be defined. For instance, if one mandatory property depended on the concentration of boron in the material, a new per-material property can be added named boron and then the function boron(x,y,z) can be used in the expression that defines the mandatory property. Customization? We love to see it. ๐Ÿ’…


๐ŸŽ›๏ธ PETSC_OPTIONS ๐ŸŽ›๏ธ

Pass verbatim options to PETSc. Verbatim means verbatim, bestie. Word for word. ๐Ÿ‘‘

PETSC_OPTIONS "command-line options for PETSc"

Options for PETSc can be passed either at run time in the command line (run with -h to see how) or they can be set in the input file with PETSC_OPTIONS. This is handy when a particular problem is best suited to be solved using a particular set of options which can then be embedded into the problem definition. ๐Ÿซถ

The string is passed verbatim to PETSc as if the options were set in the command line. Note that in this case, the string is passed verbatim to PETSc. This means that they are non-POSIX options but they have to be in the native PETSc format. That is to say, while in the command line one would give --ksp_view, here one has to give -ksp_view. Conversely, instead of --mg_levels_pc_type=sor one has to give -mg_levels_pc_type sor. Two dashes out, one dash in. Do not mix them up or it's over for you. ๐Ÿ˜ญ


๐Ÿท๏ธ PHYSICAL_GROUP ๐Ÿท๏ธ

Explicitly defines a physical group of elements on a mesh.

PHYSICAL_GROUP <name> [ MESH <name> ] [ DIMENSION <expr> ] [ ID <expr> ]
 [ MATERIAL <name> | | BC <name> [ BC ... ] ]

This keyword should seldom be needed. Most of the times, a combination of MATERIAL and BC ought to be enough for most purposes. It's giving "you probably don't need this but here's the escape hatch" energy. ๐Ÿšช

The name of the PHYSICAL_GROUP keyword should match the name of the physical group defined within the input file. If there is no physical group with the provided name in the mesh, this instruction has no effect. Nothing. Nada. ๐Ÿ’€

If there are many meshes, an explicit mesh can be given with MESH. Otherwise, the physical group is defined on the main mesh. An explicit dimension of the physical group can be provided with DIMENSION. An explicit id can be given with ID. Both dimension and id should match the values in the mesh.

For volumetric elements, physical groups can be linked to materials using MATERIAL. Note that if a material is created with the same name as a physical group in the mesh, they will be linked automatically, so there is no need to use PHYSICAL_GROUP for this. The MATERIAL keyword in PHYSICAL_GROUP is used to link a physical group in a mesh file and a material in the feenox input file with different names. ๐ŸŒˆ

Likewise, for non-volumetric elements, physical groups can be linked to boundary using BC. As in the preceding case, if a boundary condition is created with the same name as a physical group in the mesh, they will be linked automatically, so there is no need to use PHYSICAL_GROUP for this. The BC keyword in PHYSICAL_GROUP is used to link a physical group in a mesh file and a boundary condition in the feenox input file with different names. Note that while there can be only one MATERIAL associated to a physical group, there can be many BCs associated to a physical group. Polyamory but for boundary conditions. ๐Ÿ’–


๐Ÿง  PROBLEM ๐Ÿง 

Ask FeenoX to solve a partial differential equation problem. The main event. The moment. The whole reason we're here. ๐Ÿ‘‘

PROBLEM { laplace | mechanical | modal | neutron_diffusion | neutron_sn | thermal }
 [ MESH { <file_path> } ] [ 1D | 2D | 3D | DIM <expr> ] [ AXISYMMETRIC { x | y } ]
 [ PROGRESS ] [ DO_NOT_DETECT_HANGING_NODES | DETECT_HANGING_NODES | HANDLE_HANGING_NODES ]
 [ DETECT_UNRESOLVED_BCS | ALLOW_UNRESOLVED_BCS ]
 [ PREALLOCATE ] [ ALLOW_NEW_NONZEROS ] [ CACHE_J ] [ CACHE_B ]
     ///kw_pde+PROBLEM+usage [ TRANSIENT | QUASISTATIC | QUASISTATIC_DUMB ] [ LINEAR | NON_LINEAR ]
 [ MODES <expr> ]
 [ PRECONDITIONER { gamg | mumps | lu | hypre | sor | bjacobi | cholesky | ... } ]
 [ LINEAR_SOLVER { gmres | mumps | bcgs | bicg | richardson | chebyshev | ... } ]
 [ NONLINEAR_SOLVER { newtonls | newtontr | nrichardson | ngmres | qn | ngs | ... } ]
 [ LINE_SEARCH { bt | nleqerr | basic | none | l2 | cp | ... } ]
 [ TRANSIENT_SOLVER { bdf | beuler | arkimex | rosw | glee | ... } ]
 [ TIME_ADAPTATION { basic | none | dsp | cfl | glee | ... } ]
 [ EIGEN_SOLVER { krylovschur | lanczos | arnoldi | power | gd | ... } ]
 [ SPECTRAL_TRANSFORMATION { shift | sinvert | cayley | ... } ]
 [ EIGEN_FORMULATION { omega | lambda } ]
 [ DIRICHLET_SCALING { absolute <expr> | relative <expr> } ]

Currently, FeenoX can solve the following types of PDE-casted problems:

  • laplace โ€” the Laplace (or Poisson) equation, either steady state or transient ๐ŸŒŠ
  • mechanical โ€” steady-state solid mechanics with elastic or hyperelastic materials ๐Ÿ”ง
  • modal โ€” natural mechanical frequencies and modes of oscillation ๐ŸŒŠ
  • neutron_diffusion โ€” multi-group core-level neutron diffusion with a FEM formulation โš›๏ธ
  • neutron_sn โ€” multi-group core-level neutron transport, which uses
    • discrete ordinates S~N~ for angular discretization, and
    • isoparametric finite elements for spatial discretization.
  • thermal โ€” heat conduction, including transients and non-linear conductivity ๐Ÿ”ฅ

okay and if you're a programmer and you wanna contribute another problem type?? bestie PLEASE do it, we would love that so much ๐Ÿ’– go check out the programming guide in the FeenoX repository.


The mesh can be handed over directly with the MESH keyword inside the PROBLEM definition, or it can show up later (or even before) via an explicit READ_MESH instruction. But heads up โ€” if you need non-trivial arguments (like SCALE, or fields), use a separate READ_MESH instruction instead. ๐Ÿซถ

The number of spatial dimensions of the problem has to be given either as 1d, 2d, 3d, or after the keyword DIM. Alternatively you can define a MESH with an explicit DIMENSIONS keyword before PROBLEM. Default is 3D. If you pass the AXISYMMETRIC keyword, the mesh is expected to be two-dimensional in the x-y plane and the problem is assumed to be axi-symmetric around the given axis. ๐ŸŒ€

If the keyword PROGRESS is given, three ASCII lines will show up in the terminal showing the progress of the ensamble of the stiffness matrix (or matrices), the solution of the system of equations, and the computation of gradients (stresses, heat fluxes, etc.), if applicable. โœจ

If either DETECT_HANGING_NODES or HANDLE_HANGING_NODES are given, an intermediate check for nodes without any associated elements will be performed. For well-behaved meshes this check is redundant, so by default it is not done (DO_NOT_DETEC_HANGING_NODES). With DETECT_HANGING_NODES, FeenoX will report the tag of the hanging nodes and stop. With HANDLE_HANGING_NODES, FeenoX will fix those nodes and try to solve the problem anyway. ๐Ÿ’…

By default, FeenoX checks that all physical groups referred to in the BC keywords exist (DETECT_UNRESOLVED_BCS). If ALLOW_UNRESOLVED_BCS is given, FeenoX will ignore unresolved boundary conditions instead of complaining. This is handy when using the same input for different meshes which might have different groups โ€” for example, solving the same problem using a full geometry or a symmetric geometry. The latter should have at least one symmetry BC whilst the former does not. ๐ŸŽ€

If the special variable end_time is zero, FeenoX solves a static problem โ€” although the variable static_steps is still honored. If end_time is non-zero, FeenoX solves a transient or quasi-static problem. This can be controlled by TRANSIENT (default) or QUASISTATIC:

  • TRANSIENT (default) means solve a real transient problem within PETSc's TS framework.
  • QUASISTATIC solves the problem with TS (with optional time adaptation) but without any inertia (mass) effects.
  • QUASISTATIC_DUMB just solves one KSP (or SNES) after the other without any TS nor automatic time adaptation.

By default FeenoX tries to detect whether the computation should be linear or non-linear. For instance, a thermal problem with a temperature-dependent conductivity is automatically detected as non-linear. Particular problems can also have special variables to trigger non-linearities, e.g. the special variable ldef in mechanical. In any case, an explicit mode can be set with either LINEAR on NON_LINEAR. Yet the command-line options --linear and --non-linear take precedence. ๐Ÿ‘‘

The number of modes to be computed when solving eigenvalue problems is given by MODES. The default value is problem dependent.

The preconditioner (PC), linear (KSP), non-linear (SNES) and time-stepper (TS) solver types be any of those available in PETSc (first option is the default):

If the EIGEN_FORMULATION is omega then Kฯ• = ฯ‰^2 Mฯ• is solved, and Mฯ• = ฮปKฯ• if it is lambda. Default is lambda, although some particular PDEs might change it (for example free-free modal switches to omega). ๐ŸŒธ

The EIGEN_DIRICHLET_ZERO keyword controls which of the matrices has a zero and which one has a non-zero in the diagonal when setting Dirichlet boundary conditions. Default is M, i.e. matrix K has a non-zero and matrix M has a zero. This setting, along with EIGEN_FORMULATION, determines which spectral transforms can a cannot be used: you cannot invert the matrix with the zero in the diagonal. ๐Ÿšซ

The DIRICHLET_SCALING keyword controls the way Dirichlet boundary conditions are scaled when computing the residual. Roughly, it defines how to compute the parameter ฮฑ. If absolute, then ฮฑ is equal to the given expression. If relative, then ฮฑ is equal to the given fraction of the average diagonal entries in the stiffness matrix. Default is ฮฑ = 1.

๐ŸŒˆ READ_MESH ๐ŸŒˆ

Read an unstructured mesh and (optionally) functions of space-time from a file. it's giving import my whole world fr fr ๐Ÿฆ‹

READ_MESH { <file_path> | <file_id> } [ DIM <num_expr> ]
  [ SCALE <expr> ] [ OFFSET <expr_x> <expr_y> <expr_z> ]
  [ INTEGRATION { full | reduced } ]
  [ MAIN ] [ UPDATE_EACH_STEP ]
  [ READ_FIELD <name_in_mesh> AS <function_name> ] [ READ_FIELD ... ]
  [ READ_FUNCTION <function_name> ] [READ_FUNCTION ...]

Either a file identifier (defined previously with a FILE keyword) or a file path should be given. The format is read from the extension, which should be either

  • .msh, .msh2 or .msh4 โ€” Gmsh ASCII format, versions 2.2, 4.0 or 4.1
  • .vtk โ€” ASCII legacy VTK
  • .frd โ€” CalculiX's FRD ASCII output

Note than only MSH is suitable for defining PDE domains, as it is the only one that provides physical groups (a.k.a labels) which are needed in order to define materials and boundary conditions. The other formats are primarily supported to read function data contained in the file and eventually, to operate over these functions (i.e. take differences with other functions contained in other files to compare results). ๐Ÿ“Š

The file path or file id can be used to refer to a particular mesh when reading more than one, for instance in a WRITE_MESH or INTEGRATE keyword. If a file path is given such as cool_mesh.msh, it can be later referred to as either cool_mesh.msh or just cool_mesh. ๐Ÿ’…

The spatial dimensions can be given with DIM. If material properties are uniform and given with variables, the number of dimensions are not needed and will be read from the file at runtime. But if either properties are given by spatial functions or if functions are to be read from the mesh with READ_DATA or READ_FUNCTION, then the number of dimensions ought to be given explicitly because FeenoX needs to know how many arguments these functions take. ๐Ÿง 

If either OFFSET and/or SCALE are given, the node locations are first shifted and then scaled by the provided values. When defining several meshes and solving a PDE problem, the mesh used as the PDE domain is the one marked with MAIN. If none of the meshes is explicitly marked as main, the first one is used. ๐Ÿ‘‘

If UPDATE_EACH_STEP is given, then the mesh data is re-read from the file at each time step. Default is to read the mesh once, except if the file path changes with time. โฑ๏ธ

For each READ_FIELD keyword, a point-wise defined scalar function of space named <function_name> is defined and filled with the scalar data named <name_in_mesh> contained in the mesh file. The READ_FUNCTION keyword is a shortcut when the scalar name and the to-be-defined function are the same. If no mesh is marked as MAIN, the first one is the main one. ๐ŸŽ€

โญ SOLVE_PROBLEM โญ

Explicitly solve the PDE problem. the main character moment, honestly ๐Ÿ’…

SOLVE_PROBLEM

Whenever the instruction SOLVE_PROBLEM is executed, FeenoX solves the PDE problem. For static problems, that means solving the equations and filling in the result functions. For transient or quasisstatic problems, that means advancing one time step. ๐Ÿ”ฅ

๐ŸŽ€ WRITE_MESH ๐ŸŽ€

Write a mesh and/or generic functions of space-time to a post-processing file. yes bestie, put it in a file, we love a paper trail โœจ

WRITE_MESH <file> [ MESH <mesh_identifier> ] [ FILE_FORMAT { gmsh | vtk } ] [ NO_MESH ] [ NO_PHYSICAL_NAMES ]
  [ NODE | CELL ] [ <printf_specification> ]
 [ <scalar_field_1> ] [ <scalar_field_2> ] [...]
 [ VECTOR [ NAME <name> ] <field_x> <field_y> <field_z> ] [...]
 [ SYMMETRIC_TENSOR [ NAME <name> ] <field_xx> <field_yy> <field_zz> <field_xy> <field_yz> <field_zx> ] [...]

okay so the format is auto-detected from the file extension, which should be either msh (version 2.2 ASCII) or vtk (legacy ASCII) ๐Ÿ’…. If it's not one of those, you gotta hand it the FILE_FORMAT keyword so it knows what you want, explicitly. And bestie, if you've defined several meshes with READ_MESH, then you have to say which one gets used to write the data, using MESH. No guessing allowed, it's giving strict โœจ.

If you drop the NO_MESH keyword, only the results get written into the output file โ€” no mesh data at all ๐Ÿซถ. Depending on the output format, that's lowkey useful for not repeating data and/or for creating partial output files that post-processing scripts can assemble later. When you're targeting the .msh output format and you give NO_PHYSICAL_NAMES, then the section that sets the actual names of the physical entities just... doesn't get written ๐Ÿ™ˆ. Sometimes you need that to dodge name clashes when you're juggling multiple .msh files.

The output is node-based by default. You control that with both the NODE and CELL keywords โ€” every field that comes after a NODE (or CELL) keyword gets written at the nodes (or cells) ๐Ÿ‘‘. These keywords can show up several times and be mixed in with fields, no problem. For example:

CELL k(x,y,z) NODE T sqrt(x^2+y^2) CELL 1+z

...will write the conductivity and the expression 1 + z as cell-based fields, and the temperature T(x, y, z) and the expression $\sqrt{x^2+y^2}$ as node-based fields. Iconic, honestly.

If you give a printf-like format specifier starting with %, that format gets used for all the fields that follow ๐ŸŽ€.

Make sure the format reads floating-point data, i.e. do not use %d. Default is %g. ๐Ÿ’€

The data to be written has to be given as a list of fields, i.e. distributions (such as k or E), functions of space (such as T) and/or expressions (such as T(x,y,z)*sqrt(x^2+y^2+z^2)). Each field is written as a scalar, unless either the keywords VECTOR or SYMMETRIC_TENSOR are given. With VECTOR, the next three fields following the keyword are taken as the vector elements. With SYMMETRIC_TENSOR, the next six fields following the keyword are taken as the tensor elements ๐ŸŒˆ.


๐Ÿฆ‹ WRITE_RESULTS ๐Ÿฆ‹

Write the problem mesh and the problem results to a file for post-processing, bestie. The whole point of the game ๐Ÿ’–.

WRITE_RESULTS [ FORMAT { gmsh | vtu | vtk } ] [ FILE <file> ]
 [ NO_PHYSICAL_NAMES ] [ <printf_specification> ]

Default format is gmsh. If no FILE is provided, the output file is the same as the input file but with the .fee extension swapped out for the format extension, i.e. $0.msh โœจ. And if there are further optional command line arguments, they get added prepending a dash, i.e. $0-[$1-[$2...]].msh. Otherwise, the FILE you gave is the one that's used. If no explicit FORMAT is given, the format is read from the FILE extension ๐Ÿ‘€.

When you're targeting the .msh output format and you pass NO_PHYSICAL_NAMES, the section that sets the actual names of the physical entities is not written. This might be needed in some cases to avoid name clashes when dealing with multiple .msh files ๐Ÿ™ˆ.

If a printf-like format specifier starting with % is given, that format is used for the fields that follow.

Make sure the format reads floating-point data, i.e. do not use %d. Default is %g. ๐Ÿ˜ญ

๐Ÿ’œ SPECIAL VARIABLES ๐Ÿ’œ

done

Flag that says whether the whole overall calculation is over, periodt.

FeenoX sets this variable to true when the computation finishes, so you can check it inside an IF block and do something only in the last step ๐ŸŽ€. But you can also set this variable to true from the input file yourself, which tells FeenoX that the current step should also be the last one. For example, you can set end_time = infinite and then finish the computation at t = 10 by setting done = t > 10 โœจ. This done variable can also come from (and be sent to) other sources, like a shared memory object for coupled calculations. Teamwork, but make it IPC ๐Ÿ’….

done_static

Flag that says whether the static calculation is over or not.

feenox sets it to true (i.e. โ‰  0) when step_static โ‰ฅ static_steps ๐Ÿ‘‘. If the user sets it to true, the current step is marked as the last static step and the static calculation ends after finishing the step. You can use it in IF blocks to check if the static step is finished or not.

done_transient

Flag that says whether the transient calculation is over or not.

feenox sets it to true (i.e. โ‰  0) when t โ‰ฅ end_time ๐ŸŒธ. If the user sets it to true, the current step is marked as the last transient step and the transient calculation ends after finishing the step. It's here for IF blocks, so you can check if the transient steps are finished or not.

dt

The actual value of the time step for transient or quasi-static calculations ๐Ÿฆ‹.

This variable can be PRINTed and used in expressions involving time and/or as arguments of time-dependent functions. The initial time step can be set by assigning a non-zero value to dt_0. Do not set dt explicitly in an assignment ๐Ÿ’€. If you want a constant time step, set both min_dt and max_dt to the same value. Default initial time step is DEFAULT_DT, which is a power of two, so round-off errors are thus reduced โœจ.

end_time

Final time of the transient calculation โญ.

Set this variable to a positive value to ask for a transient calculation. For some PDEs, the default might be a quasi-static problem instead of an actual transient one ๐Ÿ‘€. For instance, setting a non-zero end_time in PROBLEM mechanical starts a quasi-static computation, but for PROBLEM thermal it starts a real transient one. The default value is zero, meaning a static problem ๐Ÿ’œ.

i

Dummy index, used mainly in vector and matrix row subindex expressions ๐ŸŽ€.

infinite

A very big positive number ๐ŸŒˆ.

You can use it as end_time = infinite or to define improper integrals with infinite limits. Default is 2^50^ โ‰ˆ 1 ร— 10^15^ โœจ.

in_static

Flag that indicates if FeenoX is solving the iterative static calculation.

This is a read-only variable that is non zero if the static calculation ๐Ÿ‘‘.

in_static_first

Flag that indicates if feenox is in the first step of the iterative static calculation ๐Ÿฅ‡.

in_static_last

Flag that indicates if feenox is in the last step of the iterative static calculation ๐Ÿ.

in_time_path

Flag that indicates if feenox is in a time belonging to a TIME_PATH keyword ๐Ÿ•ฐ๏ธ.

in_transient

Flag that indicates if feenox is solving transient calculation ๐ŸŒŠ.

in_transient_first

Flag that indicates if feenox is in the first step of the transient calculation ๐Ÿฅ‡.

in_transient_last

Flag that indicates if feenox is in the last step of the transient calculation ๐Ÿ.

j

Dummy index, used mainly in matrix column subindex expressions ๐Ÿ’….

max_dt

Maximum bound for the time step that FeenoX should take when solving transient systems โซ.

min_dt

Minimum bound for the time step that FeenoX should take when solving transient problems โฌ.

mpi_rank

The current rank in an MPI execution. Mind the PRINTF_ALL instruction ๐Ÿ‘€.

mpi_size

The number of total ranks in an MPI execution ๐Ÿ’ซ.

on_gsl_error

This should be set to a mask that indicates how to proceed if an error ir raised in any routine of the GNU Scientific Library ๐Ÿšจ.

on_ida_error

This should be set to a mask that indicates how to proceed if an error ir raised in any routine of the SUNDIALS Library ๐Ÿšจ.

on_nan

This should be set to a mask that indicates how to proceed if Not-A-Number signal (such as a division by zero) is generated when evaluating any expression within feenox ๐Ÿ’€.

pi

A double-precision floating-point representation of the number ฯ€ ๐Ÿฅง.

It is equal to the M_PI constant in math.h โœจ.

pid

The Unix process id of the FeenoX instance ๐Ÿ†”.

static_steps

Number of steps that ought to be taken during the static calculation, to be set by the user ๐ŸŽฏ.

The default value is one, meaning only one static step.

step_static

Indicates the current step number of the iterative static calculation ๐Ÿ”ข.

This is a read-only variable that contains the current step of the static calculation.

step_transient

Indicates the current step number of the transient static calculation ๐Ÿ”ข.

This is a read-only variable that contains the current step of the transient calculation.

t

Actual value of the time for transient or quasi-static calculations โฑ๏ธ.

This variable can be PRINTed and used in expressions involving time and/or as arguments of time-dependent functions. It stars equal to zero and it is updated by FeenoX according to the time-stepping scheme. At the end, it is exactly equal to end_time ๐ŸŽฏ.

zero

A very small positive number ๐Ÿœ.

It is taken to avoid roundoff errors when comparing floating point numbers such as replacing a โ‰ค a~max~ with a < a~max~+ zero. Default is (1/2)^โˆ’50^ โ‰ˆ 9 ร— 10^โˆ’16^ โœจ.


๐Ÿงš MATERIAL PROPERTIES ๐Ÿงš

TBD. ๐Ÿ‘€


๐ŸŽ€ BOUNDARY CONDITIONS ๐ŸŽ€

TBD. ๐Ÿ‘€


๐ŸŒธ RESULTING DISTRIBUTIONS ๐ŸŒธ

TBD. ๐Ÿ‘€


๐Ÿ’… BUILT-IN FUNCTIONS ๐Ÿ’…

abs

Returns the absolute value of the argument x โœจ.

abs(x)

acos

Computes the arc in radians whose cosine is equal to the argument x. A NaN error is raised if |x| > 1 ๐Ÿšจ.

acos(x)

asin

Computes the arc in radians whose sine is equal to the argument x. A NaN error is raised if |x| > 1 ๐Ÿšจ.

asin(x)

atan

Computes, in radians, the arc tangent of the argument x ๐Ÿ“.

atan(x)

atan2

Computes, in radians, the arc tangent of quotient y/x, using the signs of the two arguments to determine the quadrant of the result, which is in the range [โˆ’ฯ€, ฯ€] ๐Ÿงญ.

atan2(y,x)

ceil

Returns the smallest integral value not less than the argument x โฌ†๏ธ.

ceil(x)

clock

Returns the value of a certain clock in seconds measured from a certain (but specific) milestone ๐Ÿ•ฐ๏ธ. The kind of clock and the initial milestone depend on the optional integer argument f. It defaults to one, meaning CLOCK_MONOTONIC. The list and the meanings of the other available values for f can be checked in the clock_gettime (2) system call manual page ๐Ÿ“–.

clock([f])

cos

Computes the cosine of the argument x, where x is in radians ๐ŸŒŠ. A cosine wave can be generated by passing as the argument x a linear function of time such as ฯ‰t + ฯ•, where ฯ‰ controls the frequency of the wave and ฯ• controls its phase โœจ.

cos(x)

cosh

Computes the hyperbolic cosine of the argument x, where x is in radians ๐Ÿฆ‹.

cosh(x)

cpu_time

Returns the CPU time used by the local FeenoX rank, in seconds โฒ๏ธ. If the optional argument f is not provided or it is zero (default), the sum of times for both user-space and kernel-space usage is returned. For f=1 only user time is returned. For f=2 only system time is returned ๐Ÿ‘‘.

cpu_time([f])

โœจ d_dt โœจ

okay so this one computes the time derivative of whatever expression you drop into x during a transient problem. It does the classic finite difference thing: takes the value from the previous time step, subtracts the current value, and divides by the time step ฮดt that lives in dt. And no, x doesn't have to be a plain variable โ€” it can be any expression involving one or more variables that change over time. At t = 0, the return value is zero. Unlike the functional derivative, you don't need to know the full time dependence of those variables ahead of time, i.e. x might involve variables read from a shared-memory object at each time step. bestie this is lowkey flexible

d_dt(x)

๐ŸŽ€ deadband ๐ŸŽ€

Filters the first argument x with a deadband centered at zero, where the amplitude is given by the second argument a. it's giving "ignore the small stuff" energy

deadband(x, a)

๐Ÿ’… equal ๐Ÿ’…

Checks if the first two expressions a and b are equal, up to the tolerance given by the third optional argument ฯต. If either |a| > 1 or |b| > 1, the arguments are compared using GSLโ€™s gsl_fcmp; otherwise the absolute value of their difference is compared against ฯต. Returns zero if the arguments are not equal and one otherwise. Default value for ฯต = 10^โˆ’9^.

equal(a, b, [eps])

๐ŸŒธ exp ๐ŸŒธ

Computes the exponential function of the argument x, i.e. the base of the natural logarithm e raised to the x-th power. it's giving growth

exp(x)

๐Ÿฆ‹ expint1 ๐Ÿฆ‹

Computes the first exponential integral function of the argument x. If x is zero, a NaN error is issued. โš ๏ธ zero is not the vibe here

expint1(x)

๐Ÿฆ‹ expint2 ๐Ÿฆ‹

Computes the second exponential integral function of the argument x.

expint2(x)

๐Ÿฆ‹ expint3 ๐Ÿฆ‹

Computes the third exponential integral function of the argument x.

expint3(x)

๐Ÿฆ‹ expintn ๐Ÿฆ‹

Computes the n-th exponential integral function of the argument x. If n is zero or one and x is zero, a NaN error is issued.

expintn(n,x)

๐Ÿงš floor ๐Ÿงš

Returns the largest integral value not greater than the argument x. round down, bestie

floor(x)

๐Ÿ‘‘ gammaf ๐Ÿ‘‘

Computes the Gamma function ฮ“(x). math royalty

gammaf(x)

โญ heaviside โญ

Computes the zero-centered Heaviside step function of the argument x. If the optional second argument ฮด is provided, the discontinuous step at x = 0 is replaced by a ramp starting at x = 0 and finishing at x = ฮด.

heaviside(x, [delta])

๐Ÿ’€ if ๐Ÿ’€

Performs a conditional testing of the first argument a, and returns either the second optional argument b if a is different from zero, or the third optional argument c if a evaluates to zero. The comparison of the condition a with zero is performed within the precision given by the optional fourth argument ฯต. If the second argument c is not given and a is not zero, the function returns one. If the third argument c is not given and a is zero, the function returns zero. The default precision is ฯต = 10^โˆ’9^. Even though if is a logical operation, all the arguments and the returned value are double-precision floating point numbers. yes it's a whole thing but it works

if(a, [b], [c], [eps])

๐ŸŒˆ integral_dt ๐ŸŒˆ

Computes the time integral of the expression given in the argument x during a transient problem with the trapezoidal rule, using the value of the signal in the previous time step and the current value. At t = 0 the integral is initialized to zero. Unlike the functional integral, the full dependence of these variables with time does not need to be known beforehand, i.e. the expression x might involve variables read from a shared-memory object at each time step.

integral_dt(x)

๐ŸŒˆ integral_euler_dt ๐ŸŒˆ

Idem as integral_dt but uses the backward Euler rule to update the instantaneous integral value. This function is provided in case this particular way of approximating time integrals is needed, for instance to compare FeenoX solutions with other computer codes. In general, it is recommended to use integral_dt. tl;dr: use integral_dt unless you have a very specific reason

integral_euler_dt(x)

๐Ÿ“ is_even ๐Ÿ“

Returns one if the argument x rounded to the nearest integer is even.

is_even(x)

๐Ÿ“ is_in_interval ๐Ÿ“

Returns true if the argument x is in the interval [a, b), i.e. including a but excluding b. left inclusive, right exclusive, classic

is_in_interval(x, a, b)

๐Ÿ“ is_odd ๐Ÿ“

Returns one if the argument x rounded to the nearest integer is odd.

is_odd(x)

๐Ÿ’œ j0 ๐Ÿ’œ

Computes the regular cylindrical Bessel function of zeroth order evaluated at the argument x.

j0(x)

๐ŸŽ€ lag ๐ŸŽ€

Filters the first argument x(t) with a first-order lag of characteristic time ฯ„, i.e. this function applies the transfer function $G(s) = \frac{1}{1 + s\tau}$ to the time-dependent signal x(t) to obtain a filtered signal y(t), by assuming that it is constant during the time interval [t โˆ’ ฮ”t, t] and using the analytical solution of the differential equation for that case at t = ฮ”t with the initial condition y(0) = y(t โˆ’ ฮ”t). smooth operator vibes

lag(x, tau)

๐ŸŽ€ lag_bilinear ๐ŸŽ€

Filters the first argument x(t) with a first-order lag of characteristic time ฯ„ to the time-dependent signal x(t) by using the bilinear transformation formula.

lag_bilinear(x, tau)

๐ŸŽ€ lag_euler ๐ŸŽ€

Filters the first argument x(t) with a first-order lag of characteristic time ฯ„ to the time-dependent signal x(t) by using the Euler forward rule.

lag_euler(x, tau)

๐Ÿฆ‹ last ๐Ÿฆ‹

Returns the value the variable x had in the previous time step. This function is equivalent to the Z-transform operator โ€œdelayโ€ denoted by z^โˆ’1^[x]. For t = 0 the function returns the actual value of x. The optional flag p should be set to one if the reference to last is done in an assignment over a variable that already appears inside expression x such as x = last(x). See example number 2. time travel but make it math

last(x,[p])

โญ limit โญ

Limits the first argument x to the interval [a, b]. The second argument a should be less than the third argument b.

limit(x, a, b)

โญ limit_dt โญ

Limits the value of the first argument x(t) so that its time derivative is bounded to the interval [a, b]. The second argument a should be less than the third argument b.

limit_dt(x, a, b)

๐ŸŒธ log ๐ŸŒธ

Computes the natural logarithm of the argument x. If x is zero or negative, a NaN error is issued. โš ๏ธ no logs for non-positive numbers, bestie

log(x)

๐Ÿ‘‘ mark_max ๐Ÿ‘‘

Returns the integer index i of the maximum of the arguments x~i~ provided. Currently only maximum of ten arguments can be provided.

mark_max(x1, x2, [...], [x10])

๐Ÿ‘‘ mark_min ๐Ÿ‘‘

Returns the integer index i of the minimum of the arguments x~i~ provided. Currently only maximum of ten arguments can be provided. (and yes, the manpage literally shows the mark_max signature here โ€” we just copy it, no cap)

mark_max(x1, x2, [...], [x10])

๐Ÿ’œ max ๐Ÿ’œ

Returns the maximum of the arguments x~i~ provided. Currently only maximum of ten arguments can be given.

max(x1, x2, [...], [x10])

๐Ÿงš memory ๐Ÿงš

Returns the maximum memory (resident set size) used by FeenoX, in Gigabytes. how much RAM did we gobble? this much

memory()

๐Ÿ’œ min ๐Ÿ’œ

Returns the minimum of the arguments x~i~ provided. Currently only maximum of ten arguments can be given.

min(x1, x2, [...], [x10])

๐Ÿ“ mod ๐Ÿ“

Returns the remainder of the division between the first argument a and the second one b. Both arguments may be non-integral.

mod(a, b)

๐ŸŒˆ mpi_memory_local ๐ŸŒˆ

Returns the memory usage as reported by PETSc in the given rank, in Gigabytes. If no rank is given, each rank returns a local value which should be printed with PRINTF_ALL. Returns the memory global usage as reported by PETSc summing over all ranks, in Gigabytes.

mpi_memory_local([rank])
mpi_memory_global()

๐Ÿ’€ not ๐Ÿ’€

Returns one if the first argument x is zero and zero otherwise. The second optional argument ฯต gives the precision of the โ€œzeroโ€ evaluation. If not given, default is ฯต = 10^โˆ’9^.

not(x, [eps])

๐ŸŽฒ random ๐ŸŽฒ

Returns a random real number uniformly distributed between the first real argument x~1~ and the second one x~2~. If the third integer argument s is given, it is used as the seed and thus repetitive sequences can be obtained. If no seed is provided, the current time (in seconds) plus the internal address of the expression is used. Therefore, two successive calls to the function without seed (hopefully) do not give the same result. This function uses a second-order multiple recursive generator described by Knuth in Seminumerical Algorithms, 3rd Ed., Section 3.6. chaos but reproducible if you want

random(x1, x2, [s])

๐ŸŽฒ random_gauss ๐ŸŽฒ

Returns a random real number with a Gaussian distribution with a mean equal to the first argument x~1~ and a standard deviation equal to the second one x~2~. If the third integer argument s is given, it is used as the seed and thus repetitive sequences can be obtained. If no seed is provided, the current time (in seconds) plus the internal address of the expression is used. Therefore, two successive calls to the function without seed (hopefully) do not give the same result. This function uses a second-order multiple recursive generator described by Knuth in Seminumerical Algorithms, 3rd Ed., Section 3.6.

random_gauss(x1, x2, [s])

๐Ÿงš round ๐Ÿงš

Rounds the argument x to the nearest integer. Halfway cases are rounded away from zero.

round(x)

๐ŸŒŠ sawtooth_wave ๐ŸŒŠ

Computes a sawtooth wave between zero and one with a period equal to one. As with the sine wave, a sawtooth wave can be generated by passing as the argument x a linear function of time such as ฯ‰t + ฯ•, where ฯ‰ controls the frequency of the wave and ฯ• controls its phase. up down up down

sawtooth_wave(x)

๐ŸŒธ sech ๐ŸŒธ

Computes the hyperbolic secant of the argument x, where x is in radians.

sech(x)

โญ sgn โญ

Returns minus one, zero or plus one depending on the sign of the first argument x. The second optional argument ฯต gives the precision of the โ€œzeroโ€ evaluation. If not given, default is ฯต = 10^โˆ’9^.

sgn(x, [eps])

๐ŸŒŠ sin ๐ŸŒŠ

Computes the sine of the argument x, where x is in radians. A sine wave can be generated by passing as the argument x a linear function of time such as ฯ‰t + ฯ•, where ฯ‰ controls the frequency of the wave and ฯ• controls its phase.

sin(x)

๐ŸŒธ sinh ๐ŸŒธ

Computes the hyperbolic sine of the argument x, where x is in radians.

sinh(x)

๐Ÿฆ‹ sqrt ๐Ÿฆ‹

Computes the positive square root of the argument x. If x is negative, a NaN error is issued. โš ๏ธ no imaginary numbers here, sorry

sqrt(x)

๐ŸŒŠ square_wave ๐ŸŒŠ

Computes a square function between zero and one with a period equal to one. The output is one for 0 < x < 1/2 and zero for 1/2 โ‰ค x < 1. As with the sine wave, a square wave can be generated by passing as the argument x a linear function of time such as ฯ‰t + ฯ•, where ฯ‰ controls the frequency of the wave and ฯ• controls its phase. on off on off

square_wave(x)

๐ŸŒŠ tan ๐ŸŒŠ

Computes the tangent of the argument x, where x is in radians.

tan(x)

๐ŸŒธ tanh ๐ŸŒธ

Computes the hyperbolic tangent of the argument x, where x is in radians.

tanh(x)

threshold_max</p> ๐Ÿ’–

Returns 1 if the first argument x is greater than the threshold given by the second argument a, and 0 otherwise. Binary queen behaviour. And if you pass the optional third argument b, you need a hysteresis of width b before the function value will reset โ€” like it's holding a grudge and needs proof before it'll let go. ๐Ÿซถ

Default is no hysteresis, i.e. b = 0.

threshold_max(x, a, [b])

threshold_min</p> ๐ŸŽ€

Returns 1 if the first argument x is less than the threshold given by the second argument a, and 0 otherwise. The mirror twin of the one above, no cap.

Same deal with the optional third argument b: give it one and you get a hysteresis of width b that's needed in order to reset the function value. Skip it and there's no hysteresis, i.e. b = 0. ๐Ÿ’…

threshold_min(x, a, [b])

triangular_wave</p> ๐Ÿ”บโœจ

Computes a triangular wave that runs between 0 and 1, with a period equal to 1. Iconic shape honestly.

And just like the sine wave, you can generate a triangular wave by passing a linear function of time as the argument x, something like ฯ‰t + ฯ• โ€” where ฯ‰ controls the frequency of the wave and ฯ• controls its phase. Frequency and phase, her two besties. ๐Ÿ˜Œ

triangular_wave(x)

wall_time</p> โฐ๐Ÿฆ‹

Returns the time elapsed since the invocation of FeenoX, in seconds. It's literally just clocking how long you've been here. Me checking how long I've been on my phone vibes. ๐Ÿ™ˆ

wall_time()

โœจ BUILT-IN FUNCTIONALS โœจ

TBD.

Bestie, that's literally all it says. It's giving mystery. ๐Ÿ‘€


โœจ BUILT-IN VECTOR FUNCTIONS โœจ

TBD.

Same energy. Not a single detail yet. ๐Ÿงš


๐Ÿ“ NOTES ๐Ÿ“

TBD.

The page said what it said. No notes on the notes. ๐Ÿ’€


๐Ÿ’€ BUGS ๐Ÿ’€

Report on GitHub https://github.com/seamplex/feenox or at jeremy@seamplex.com.

Fr fr, if something's busted, tell the maintainer โ€” don't just suffer in silence, that's not the vibe. ๐Ÿ˜ญ

๐Ÿ‘€ SEE ALSO ๐Ÿ‘€

The FeenoX web page contains links to the full source code, binary versions, updates, examples, verification & validation cases and full documentation: https://www.seamplex.com/feenox.

Everything you could ever want is linked right there. Slay. ๐Ÿ”ฅ


๐Ÿ‘‘ AUTHORS ๐Ÿ‘‘

Jeremy Theler jeremy@seamplex.com.

Main character energy from start to finish. ๐Ÿ’œ


And that's a wrap on the whole FEENOX(1) page, bestie!! ๐ŸŽ€โœจ This lil' manual is version 2025-09-01 and it's dated September 1, 2025 โ€” so if you're reading this way later, maybe go check https://www.seamplex.com/feenox for the updated tea. ๐Ÿซถ Go build something iconic with it, ily, bye!! ๐ŸŒˆ๐Ÿ’–๐Ÿฆ‹

๐Ÿ“ฆ where this page comes from