Started work on mkdocs documentation site.

docs/febioheat_theory.md

@@ -0,0 +1,126 @@
+## Introduction
+Theory Manual
+The heat transfer equation is defined as follows.
+
+\[
+\begin{equation}
+\label{eq:heat-transfer}
+\rho C\frac{\partial u}{\partial t}+\nabla \cdot \mathbf{q}=Q
+\end{equation}
+\]
+
+Here, $u$ is the _temperature_, $\rho$ the _density_, $C$ the _specific heat_, $\mathbf{q}$ the _heat flux_, and $Q$ the _heat source_. 
+
+The following boundary conditions are assumed.
+
+\[\begin{align}
+  & u={{u}_{0}}\quad\text{on}\,{{\Gamma }_{D}} \\ 
+ & \mathbf{q}\cdot \mathbf{n}={{q}_{0}}\quad\text{on }\,{{\Gamma }_{N}} \\ 
+\end{align}\]
+
+
+Furthermore, it is assumed that the heat flux is defined via Fourier’s law,
+
+\[
+\begin{equation}
+\label{eq:heat-flux}
+\mathbf{q}=-k\,\nabla u
+\end{equation}
+\]
+
+where $k$ is the thermal conductivity. 
+
+## Weak Formulation
+To construct the weak form, first multiply the left- and right-hand side with an arbitrary function $w$ (that satisfies the homogenous boundary conditions, i.e. $w=0$ on ${{\Gamma }_{D}}$).
+
+\[w\left( \rho C\frac{\partial u}{\partial t}+\nabla \cdot \mathbf{q} \right)=wQ\]
+
+Then, integrate over the problem domain.
+
+\[\int\limits_{\Omega }{w\left( \rho C\frac{\partial u}{\partial t}+\nabla \cdot \mathbf{q} \right)d\Omega }=\int\limits_{\Omega }{wQd\Omega }\]
+
+Using integration by parts and the boundary conditions, we can rewrite the second term on the left-hand side.
+
+\[\int\limits_{\Omega }{w\rho C\frac{\partial u}{\partial t}d\Omega }+\int\limits_{{{\Gamma }_{N}}}{w{{q}_{0}}d\Gamma }-\int\limits_{\Omega }{\nabla w\cdot \mathbf{q}d\Omega }=\int\limits_{\Omega }{wQd\Omega }\]
+
+Using equation (\ref{eq:heat-flux}) results in,
+
+\[\int\limits_{\Omega }{w\rho C\frac{\partial u}{\partial t}d\Omega }+\int\limits_{\Omega }{k\nabla w\cdot \nabla ud\Omega }=\int\limits_{\Omega }{wQd\Omega }-\int\limits_{{{\Gamma }_{N}}}{w{{q}_{0}}d\Gamma }\]
+
+For steady-state analysis, we simply omit the first term.
+
+\[\int\limits_{\Omega }{k\nabla w\cdot \nabla ud\Omega }=\int\limits_{\Omega }{wQd\Omega }-\int\limits_{{{\Gamma }_{N}}}{w{{q}_{0}}d\Gamma }\quad\text{(steady state analysis)}\]
+
+## Discretization
+The domain is discretized using iso-parametric elements.
+
+\[\begin{align}
+  & \mathbf{r}=\sum\limits_{a}{{{N}_{a}}{{\mathbf{r}}_{a}}} \\ 
+ & u=\sum\limits_{a}{{{N}_{a}}{{u}_{a}}} \\ 
+ & w=\sum\limits_{a}{{{N}_{a}}{{w}_{a}}} \\ 
+\end{align}\]
+
+Considering first the simpler case of steady state, the discretized equation becomes the following.
+
+\[\sum\limits_{a,b}{{{w}_{a}}\left( \int\limits_{\Omega }{k\nabla {{N}_{a}}\cdot \nabla {{N}_{b}}d\Omega } \right){{u}_{b}}}=\sum\limits_{a}{{{w}_{a}}\left( \int\limits_{\Omega }{{{N}_{a}}Qd\Omega } \right)}-\sum\limits_{a}{{{w}_{a}}\left( \int\limits_{{{\Gamma }_{N}}}{N_{a}^{'}{{q}_{0}}d\Gamma } \right)}\]
+
+Note that in the second term on the right-hand side, the shape functions are restricted to the boundary surface.
+
+This can be written more concisely.
+
+\[\mathbf{Ku}=\mathbf{F}\]
+
+Here, $\mathbf{K}$ is the stiffness matrix and is calculated from,
+
+\[{{K}_{ab}}=\int\limits_{\Omega }{k\nabla {{N}_{a}}\cdot \nabla {{N}_{b}}d\Omega }\]
+
+Note that the stiffness matrix is symmetric.
+
+The vector $\mathbf{u}$ collects all the unknown temperature values, and $\mathbf{F}$ collects all the known force loads,
+
+\[{{F}_{a}}=\int\limits_{\Omega }{{{N}_{a}}Qd\Omega }-\int\limits_{{{\Gamma }_{N}}}{N_{a}^{'}{{q}_{0}}d\Gamma }\]
+
+Note again the slight abuse of notation, since the second term is evaluated over a surface, using different shape functions than the first term.
+
+Returning now to the transient case, we still need to discretize the term,
+
+\[\int\limits_{\Omega }{w\rho C\dot{u}d\Omega }\]
+
+Here, we use $\dot{u}$ to denote the time derivative of the temperature. Introducing shape functions, similar as before,
+
+\[\dot{u}=\sum\limits_{a}{{{N}_{a}}{{{\dot{u}}}_{a}}}\]
+
+The discretized integral then becomes,
+
+\[\sum\limits_{a,b}{{{w}_{a}}\left( \int\limits_{\Omega }{\rho C{{N}_{a}}{{N}_{b}}d\Omega } \right)\,{{{\dot{u}}}_{b}}}=\sum\limits_{a,b}{{{w}_{a}}{{M}_{ab}}{{{\dot{u}}}_{b}}}\]
+
+where,
+
+\[{{M}_{ab}}=\int\limits_{\Omega }{\rho C{{N}_{a}}{{N}_{b}}d\Omega }\,\]
+
+is the _capacitance matrix_.
+
+Bringing all the terms together results in the following equation,
+
+\[\mathbf{M\dot{u}}+\mathbf{Ku}=\mathbf{F}\]
+
+This is the _semi-discrete equation_. Solving it will be discussed below.
+
+
+## Time Discretization
+
+We now turn our attention to solving the semi-discrete equation, which is repeated here.
+
+\[\mathbf{M\dot{u}}+\mathbf{Ku}=\mathbf{F}\]
+
+We’ll solve this equation using the implicit-Euler method (also known as backward Euler method). This method is unconditionally stable and first-order accurate.
+
+First, we replace $\dot{\mathbf{u}}$ with a finite difference approximation.
+
+\[\mathbf{\dot{u}}=\frac{1}{\Delta t}\left( {{\mathbf{u}}^{n+1}}-{{\mathbf{u}}^{n}} \right)\]
+
+Then, we solve the semi-discrete equation at time point n+1.
+
+\[\left( \frac{1}{\Delta t}\mathbf{M}+\mathbf{K} \right){{\mathbf{u}}^{n+1}}=\mathbf{F}+\frac{1}{\Delta t}\mathbf{M}{{\mathbf{u}}^{n}}\]
+
+Once again, we recover a linear system of equations that can be solved.

docs/febioheat_tutorial1.md

@@ -0,0 +1,3 @@
+# Tutorial 1
+
+(coming soon to a manual nearby!)

docs/febioheat_user_manual.md

@@ -0,0 +1,128 @@
+## Introduction
+
+This user manual discusses the feature of the FEBioHeat plugin in the context of the FEBio input file. 
+
+This manual assumes the 4.0 FEBio format specification.
+
+## Heat Module
+
+The module type must be set to `heat` for solving a heat-transfer problem with this plugin. 
+
+```xml
+<Module type="heat"/>
+```
+
+## Heat transfer constitutive model
+
+Currently, only one constitutive model is supported for heat-transfer analysis, namely the `isotropic Fourier`, which defines the heat flux $\mathbf{q}$ as follows,
+
+\[\mathbf{q}=-k\,\nabla T\]
+
+Here, $T$ is the temperature and $k$ the _thermal conductivity_.
+
+This material has the following parameters as shown in the following table. 
+
+| parameter | description               |units (SI)|
+|-----------|---------------------------|----------|
+|`density`  | The material density      | kg/m^3   |
+|`k`        | The thermal conductivity	| W/m.K    |
+|`c`        | The specific heat         | J/kg.K   |
+
+The units are just given as an example of a consistent set of units for these variables. Users can use their preferred unit system instead. 
+
+Note that density and specific heat are only used for transient analysis.
+
+An example of material definition follows below.
+
+```xml
+<material id="1" name="myMaterial1" type="isotropic Fourier">
+	<density>1.0</density>
+	<k>0.4</k>
+	<c>1.0</c>
+</material>
+```
+
+## Heat Transfer Boundary Conditions
+### Fixed and prescribed temperature
+
+To prescribe a temperature on a boundary, the “prescribed temperature” boundary condition can be used. A load controller can be assigned to the value parameter to make the prescribed temperature a function of time. For example,
+
+```xml
+<bc node_set="PrescribedBC1" type="prescribed temperature">
+	<value lc="1">1.0</value>
+	<relative>0</relative>
+</bc>
+```
+
+For prescribing a zero temperature on a boundary, a slightly more efficient boundary condition can be used instead, namely the `zero temperature` boundary condition.
+
+```xml
+<bc name="BC01" node_set="FixedTemperature1" type="zero temperature"/>
+```
+
+This is more efficient since the degree of freedom for nodes assigned to a `zero temperature` boundary condition, is effectively ignored, resulting in smaller linear system of equations that needs to be solved.
+
+### Heat flux
+The heat flux surface load prescribes the heat flux on a surface of the mesh. It is defined as a surface_load with the type attribute set to `heatflux`. It takes a single parameter, namely `flux`, which specifies the value of the heat flux.
+
+```xml
+<surface_load type="heatflux" surface="Surface01">
+	<flux lc="1">2.5</flux>
+</surface_load>
+```
+
+### Convective heat flux
+A convective heat flux applies a flux boundary condition where the flux is proportional to the difference between the surface temperature $T$ and the ambient temperature $T_a$.
+
+\[{{q}_{c}}={{h}_{c}}\left( T-{{T}_{a}} \right)\]
+
+This boundary load is defined as a `surface_load` with the type attribute set to `convective_heatflux` and requires two parameters.
+
+| parameter | description                  | units (SI)|
+|-----------|------------------------------|-----------|
+| `hc`      | The proportionality constant | W/m2.K    |
+| `Ta`      | The ambient temperature      | K         |
+
+An example is given below. 
+
+```xml
+<surface_load type="convective_heatflux" surface="Surface01">
+	<hc>60.0</hc>
+	<Ta lc="1">25</Ta>
+</surface_load>
+```
+
+## Heat source
+
+A heat source can be added by defining a body_load with the type set to `heat_source`. Only one parameter is required, `Q`, to define the heat source value.
+
+| parameter | description                    | units (SI)|
+|-----------|--------------------------------|-----------|
+| `Q`       | The constant heat source value | W/m^3      |
+
+```xml
+<body_load type="heat_source">
+	<Q lc="1">13.5</Q>
+</body_load>
+```
+
+## Output variables
+
+The heat transfer plugin adds several output variables to both the plot file and the log file.
+
+### Plot file variables
+The following table shows a list of all new variables defined by the plugin.
+
+| variable    | description                           |
+|-------------|---------------------------------------|
+|`temperature`|The nodal temperature at the mesh nodes|
+|`heat flux`  |The average heat flux over each element|
+
+### Log file variables
+
+The following table shows a list a log file variables.
+
+|variable | description       |class    |
+|---------|-------------------|---------|
+| `T`     | nodal temperature |node_data|
+

docs/index.md

@@ -0,0 +1,14 @@
+---
+date:
+  created: 2026-1-18
+---
+
+# FEBioHeat
+
+This is the documentation for **FEBioHeat**, an FEBio plugin that allows users to solve steady-state and transient heat transfer problems. 
+
+Please see the [Tutorial](febioheat_tutorial1.md) page for a practical guide on using the FEBioHeat plugin with FEBioStudio. 
+
+The [User Manual](febioheat_user_manual.md) page describes the features of the plugin in the context of the FEBio Input file. 
+
+For some background on the theory, please see the [Theory](febioheat_theory.md) page. 

docs/js/mathjax_config.js

@@ -0,0 +1,7 @@
+window.MathJax = {
+  tex: {
+    tags: 'ams',
+    inlineMath: [['$', '$'], ['\\(', '\\)']],
+    displayMath: [['$$', '$$'], ['\\[', '\\]']]
+  }
+};

mkdocs.yml

@@ -0,0 +1,51 @@
+site_name: FEBioHeat
+site_description: Documentation site for FEBioHeat
+site_author: Steve Maas
+theme:
+  name: material
+  logo: figs/FEBioHeat.png
+  palette:
+    primary: indigo
+    accent: indigo
+  font:
+    text: 'Roboto'
+    code: 'Roboto Mono'
+  features:
+    - navigation.tabs
+    - navigation.top
+    - search.highlight
+    - search.suggest
+    - toc.integrate
+    - content.external.links
+    - content.code.annotate
+    - content.code.copy
+markdown_extensions:
+  - admonition
+  - codehilite
+  - footnotes
+  - attr_list
+  - md_in_html  
+  - toc:
+      permalink: true
+  - pymdownx.arithmatex:
+      generic: true
+  - pymdownx.details
+  - pymdownx.superfences
+  - pymdownx.blocks.caption
+  - def_list
+  - pymdownx.tasklist:
+      custom_checkbox: true
+  - pymdownx.emoji:
+      emoji_index: !!python/name:material.extensions.emoji.twemoji
+      emoji_generator: !!python/name:material.extensions.emoji.to_svg
+
+extra_javascript:
+  - js/mathjax_config.js
+  - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js
+
+nav:
+  - Home: index.md
+  - Tutorials:
+      - Tutorial 1: febioheat_tutorial1.md
+  - User Manual: febioheat_user_manual.md
+  - Theory: febioheat_theory.md