diff --git a/content/project/femformal/featured.png b/content/project/femformal/featured.png
new file mode 100644
index 0000000..f5a153a
Binary files /dev/null and b/content/project/femformal/featured.png differ
diff --git a/content/project/femformal/index.md b/content/project/femformal/index.md
new file mode 100644
index 0000000..29a235a
--- /dev/null
+++ b/content/project/femformal/index.md
@@ -0,0 +1,206 @@
+---
+title: Formal Methods for PDEs
+summary: A Formal Methods Approach to Synthesis and Verification of Systems Governed by Partial Differential Equations
+tags: []
+date: "2021-02-18"
+
+# Optional external URL for project (replaces project detail page).
+external_link: ""
+
+image:
+  caption: ""
+  focal_point: Smart
+
+links:
+# - icon: twitter
+#   icon_pack: fab
+#   name: Follow
+#   url: https://twitter.com/georgecushen
+url_code: "https://github.com/fran-penedo/femformal"
+url_pdf: ""
+url_slides: ""
+url_video: ""
+
+# Slides (optional).
+#   Associate this project with Markdown slides.
+#   Simply enter your slide deck's filename without extension.
+#   E.g. `slides = "example-slides"` references `content/slides/example-slides.md`.
+#   Otherwise, set `slides = ""`.
+slides: ""
+---
+Partial differential equations (PDEs) model nearly all of the physical systems
+and processes of interest to scientists and engineers. Some well-known examples
+include the Navier-Stokes equation for fluid mechanics, the Maxwell equations
+for electromagnetics, the Schrödinger equation for quantum mechanics, the heat
+equation and the wave equation. The analysis of these and other PDEs has had a
+tremendous impact on society by enabling our understanding of thermal,
+electrical, fluidic and mechanical processes.
+
+While a mature field, the study of PDEs is often approached through simulations
+in which approximate models obtained through spatio-temporal
+discretization techniques, such as the Finite Element Method (FEM),
+are solved numerically. These methods, however,
+do not allow for rigorous guarantees that a system satisfies a complex
+specification. Other approaches generalizing classic Ordinary Differential
+Equations (ODEs) tools such as Lyapunov
+analysis and backstepping do provide some guarantees
+and can be used to obtain boundary control strategies for a wide variety of
+PDEs. However, they are restricted to classical control objectives such as
+stabilization and cost optimization.
+
+The formal statement of specifications and the development of
+analysis techniques that can guarantee their satisfaction by design has been the main
+focus of the formal methods field. During the past decades, many specification
+languages have been defined, such as Linear Temporal
+Logic (LTL) and Signal Temporal Logic (STL).
+Originally, these logics have been applied to the
+study of finite systems, although more recently
+abstraction procedures have been developed to reduce problems with ODEs
+to finite models (for example through state space
+discretization or mixed-integer linear program (MILP) encodings).
+However, these techniques cannot be immediately applied to the analysis of PDEs
+due to the lack of spatial information in the formal language. This issue was
+first addressed in in the context of spatially
+distributed systems, which can be viewed as PDEs in a discretized spatial domain.
+In their work, the authors view the system state as an image and define a formal
+language with explicit spatial information called Linear Superposition Logic (LSSL)
+based on quadtrees, a tree-based
+abstraction of the system. A variant of LSSL with quantitative semantics, 
+Tree Spatial Superposition Logic (TSSL), was used 
+for (steady-state) pattern synthesis in reaction
+diffusion systems, and a formal
+language combining TSSL and STL, Spatial Temporal Logic (SpaTeL), allows the
+synthesis of dynamical patterns. The specification of patterns in these logics
+is difficult, however, and machine learning techniques are needed in order to
+obtain logic formulas corresponding to a set of examples of the desired pattern,
+which is not always available or desirable for some applications where
+the system behavior must be specified a priori. More recently, a new
+spatio-temporal logic called STREL was introduced in
+to specify the behavior of mobile and
+spatially distributed Cyber-Physical Systems.
+
+One of the major drawbacks of these methods is the issue of scalability. Both of the
+main approaches, automata and optimization, quickly become unfeasible when the system dimension
+increases, which poses a problem when dealing with infinite-dimensional systems such as
+PDEs. In the case of optimization based techniques, a recent trend has focused on
+decentralized synthesis and verification where the coupling between subsystems is
+formally captured in Assume-Guarantee Contracts (AGCs). However, the key
+assumption of monotonicity is not present in PDE systems.
+
+Besides classical problems that focus on the
+temporal evolution of the state of the system, PDEs allow for a different kind of
+problems where the aim is to specify the intrinsic behavior of the system through its so
+called constitutive response. In this case, instead of designing boundary control
+inputs, the engineer is tasked with the design of the geometry of a structure that
+exhibits the desired mechanical or thermal properties. The traditional
+approach to this problem, topology optimization is not enough to tackle
+both geometric and material nonlinearity.
+
+## Boundary Control
+
+In order to solve the classical boundary control problems through the use of formal
+methods, we've developed a formal language called Spatial-STL (S-STL) and a framework
+that reformulates a verification or control synthesis problem for a PDE system into an
+optimization problem for a system of difference equations that can be encoded and solved
+as a Mixed-Integer Linear Program (MILP). As an example of the kind of properties we
+can describe with S-STL, consider the following: 
+
+Consider a metallic rod of 100mm. The
+temperature at one end of the rod is fixed at 300 kelvin, while a
+heat source is applied to the other end. The temperature of the rod follows
+a linear 1D heat equation. We want the temperature distribution of the
+rod to be within 3 kelvin of the linear profile $\mu(x) = \frac{x}{4} + 300$ 
+at all times between 4 and 5 seconds in the section between 30 and 60
+mm. Furthermore, the temperature should never exceed 345 kelvin
+at the point where the heat source is applied ($x=100$).
+We can formulate such a specification using the following S-STL formula:
+
+\begin{equation}\label{eq:phi_ex}
+    \begin{aligned}
+    \phi_{ex} = &\mathbf{G}_{[4,5]} 
+    \big((\forall x \in [30,60] : u(x) - (\frac{x}{4} + 303) < 0 )  \land \\\\
+          &\quad (\forall x \in [30, 60] : u(x) - (\frac{x}{4} + 297) >
+        0)\big) \land \\\\
+      &\mathbf{G}_{[0, 5]} (\forall x \in [100, 100] : u(x) -
+    345 < 0)
+    \end{aligned}
+\end{equation}
+
+Note that S-STL admits quantitative semantics similar to STL (see the [templogic project
+page](/project/templogic) for a quick introduction.)
+
+We implemented this framework in a Python tool
+called [FEMFormal](https://github.com/fran-penedo/femformal). In FEMFormal, you can
+define boundary control problems in 1D and 2D linear heat and wave equations, as well as
+some non-linear 1D wave equations. Then, you can solve verification and synthesis
+problems with S-STL specifications such as the one above. For example, for a mixed steel
+and brass rod, our tool obtains the following solution:
+
+![Heat equation synthesis trajectory](/media/femformal_heat.gif)
+![Heat equation synthesis inputs](/media/femformal_heat_inputs.png)
+
+You can see more solved examples as well as how to use FEMFormal
+[here](https://github.com/fran-penedo/femformal).
+
+## Constitutive Response Design
+
+Suppose we now want to produce a wide variety of mechanical behaviors out of a
+single material. For example, can we give an elastic material a tunable,
+effective yield point, or make it strain harden or soften at a specified amount
+of displacement? This response to loading can be determined from the
+stress-strain curve, which serves as a fundamental piece of information that
+describes the global mechanical properties of structures and materials. Recent
+advances in mechanical metamaterials have shown that it is possible to tailor
+the relationship between stress and strain of a particular material through
+structural patterning.
+
+{{< figure src="/media/femformal_bertoldi.png" title="Metamaterial design (from Shim et al., 2013)" >}}
+
+In order to formally specify the desired response
+of the structure, we developed a formal language over stress-strain curves and
+adapted S-STL for this purpose. Then, we proposed an optimization
+procedure based on gradient-free optimization algorithms and simulation.
+
+Consider the following example of this process: a square sheet of rubber
+of side $L = 8cm$ and thickness $h = 1cm$ is drilled with
+circular holes of radius $r$ constrained to be within the interval
+$[0.3cm, 0.4cm]$ and distributed in a lattice given by the
+vectors $v_1, v_2 \in \mathbb{R}^2$ such that they do not intersect. 
+The constitutive response was obtained as a
+positive force-displacement curve by loading the material with a force $F(t)$ applied
+uniformly across the top surface such that the discplacement at the top surface
+increases linearly from 0 to 1.05cm in 3.5 seconds. 
+The specification in this case is to have the
+force-displacement curve of the material to be within $\epsilon$ tolerance of a
+target curve $p$, which can be expressed in S-STL as the following formula:
+
+\begin{equation}
+\begin{aligned}
+    \phi = \forall x \in \{(0, L)\} : \mathbf{G}_{[0, T]} 
+        (F < p(d_2(x)) + \epsilon \land F > p(d_2(x)) - \epsilon) \\\\
+        p(d) = \begin{cases}
+            \frac{40}{7.5 \cdot 10^{-3}} d & d < 7.5 \cdot 10^{-3} \\\\
+            40 & d > 7.5 \cdot 10^{-3} 
+         \end{cases}
+\end{aligned}
+\end{equation}
+
+where $d_2(x)$ represents the displacement of the top boundary at $x$, $F$ is
+the applied force, the temporal operator $\mathbf{G}_{[0, T]} \psi$ means that
+$\psi$ must be satisfied at all times in the interval and the spatial 
+operator $\forall x \in \{(0, L)\} : \psi$ means that $\psi$ must be satisfied at all points in the
+interval.
+
+We solved this problem using differential evolution to find the best pattern
+with respect to the S-STL score of the specification
+$\phi$. We generate simulations of the system using the Finite Element Method simulator 
+[ANSYS](https://www.ansys.com/).
+After 280 evaluations performed in about 38 hours, we obtained the following pattern:
+
+![Constitutive Response Results](/media/const_design1.png)
+![Constitutive Response Results](/media/const_design2.png)
+
+This pattern results in a S-STL score of 0.641, with the following force-displacement
+curve (compared to the target curve and its tolerance):
+
+![Constitutive Response F-D Curve](/media/const_design_response.png)
diff --git a/content/publication/cdc2018/index.md b/content/publication/cdc2018/index.md
index 079c98c..56417ca 100644
--- a/content/publication/cdc2018/index.md
+++ b/content/publication/cdc2018/index.md
@@ -67,7 +67,7 @@ image:
 #   E.g. `internal-project` references `content/project/internal-project/index.md`.
 #   Otherwise, set `projects: []`.
 projects:
-- formal-pde
+- femformal
 - templogic
 
 # Slides (optional).
diff --git a/content/publication/thesis/featured.png b/content/publication/thesis/featured.png
index e5458e0..f5a153a 100644
Binary files a/content/publication/thesis/featured.png and b/content/publication/thesis/featured.png differ
diff --git a/content/publication/thesis/index.md b/content/publication/thesis/index.md
index 2ae87bb..cca7994 100644
--- a/content/publication/thesis/index.md
+++ b/content/publication/thesis/index.md
@@ -74,7 +74,7 @@ image:
 #   E.g. `internal-project` references `content/project/internal-project/index.md`.
 #   Otherwise, set `projects: []`.
 projects:
-- formal-pde
+- femformal
 - templogic
 - agmilp
 
diff --git a/static/media/const_design1.png b/static/media/const_design1.png
new file mode 100644
index 0000000..ed223b9
Binary files /dev/null and b/static/media/const_design1.png differ
diff --git a/static/media/const_design2.png b/static/media/const_design2.png
new file mode 100644
index 0000000..c00fffb
Binary files /dev/null and b/static/media/const_design2.png differ
diff --git a/static/media/const_design_response.png b/static/media/const_design_response.png
new file mode 100644
index 0000000..158bf60
Binary files /dev/null and b/static/media/const_design_response.png differ
diff --git a/static/media/femformal.gif b/static/media/femformal.gif
new file mode 100644
index 0000000..22b1b07
Binary files /dev/null and b/static/media/femformal.gif differ
diff --git a/static/media/femformal_bertoldi.png b/static/media/femformal_bertoldi.png
new file mode 100644
index 0000000..bd20351
Binary files /dev/null and b/static/media/femformal_bertoldi.png differ
diff --git a/static/media/femformal_heat.gif b/static/media/femformal_heat.gif
new file mode 100644
index 0000000..052ee33
Binary files /dev/null and b/static/media/femformal_heat.gif differ
diff --git a/static/media/femformal_heat_inputs.png b/static/media/femformal_heat_inputs.png
new file mode 100644
index 0000000..6330bbd
Binary files /dev/null and b/static/media/femformal_heat_inputs.png differ