/*

Copyright (c) 2005-2024, University of Oxford.

All rights reserved.

University of Oxford means the Chancellor, Masters and Scholars of the

University of Oxford, having an administrative office at Wellington

Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice,

this list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,

this list of conditions and the following disclaimer in the documentation

and/or other materials provided with the distribution.

* Neither the name of the University of Oxford nor the names of its

contributors may be used to endorse or promote products derived from this

software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE

LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR

CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE

GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)

HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT

LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT

OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

/*

*

*

*

* Chaste tutorial - this page gets automatically changed to a wiki page

* DO NOT remove the comments below, and if the code has to be changed in

* order to run, please check the comments are still accurate

*

*

*

*

*

*/

#ifndef TESTWRITINGPDESOLVERSTUTORIAL_HPP_

#define TESTWRITINGPDESOLVERSTUTORIAL_HPP_

/* HOW_TO_TAG PDE

* Write new PDE solvers (especially for linear coupled elliptic/parabolic systems)

*/

/*

* ### Introduction

*

* Chaste can be used to set up solvers for more general (coupled) PDEs. To do this the

* user just needs to code up the integrands of any finite element (FE) matrices or vectors,

* without having to deal with set-up, looping over elements, numerical quadrature, assembly

* or solving the linear system. If you have a set of coupled linear PDEs for which it is

* appropriate to use linear basis functions for each unknown (for example, a reaction-diffusion

* system), then it is relatively straightforward to set up a solver that will be parallel and

* reliable, since all the base components are heavily tested.

*

* Some solvers for general simple (uncoupled) linear PDEs are already provided in Chaste, such

* as the `SimpleLinearEllipticSolver`. These are for PDEs that can be written in a generic

* form (`SimpleLinearEllipticPde`, for example). However ''coupled'' PDEs can't

* be easily written in generic form, so the user has to write their own solver. In this tutorial

* we explain how to do this.

*

* For this tutorial the user needs to have read the solving-PDEs tutorials. It may also be

* helpful to read the associated [lectures notes](https://chaste.cs.ox.ac.uk/trac/wiki/ChasteGuides/NmoodLectureNotes),

* in particular the slides on solving equations using finite elements if you are not familiar

* with this (lecture 2), the slides on the general design of the Chaste finite element solvers

* (lecture 3), and the first part of lecture 4.

*

* Let us use the terminology "assembled in an FE manner" for any matrix or vector that is

* defined via a volume/surface/line integral, and which is constructed by: looping over

* elements (or surface elements, etc); computing the elemental contribution (i.e. a small

* matrix/vector) using numerical quadrature; and adding to the full matrix/vector.

*

* We only consider linear problems here. In these problems the discretised FE problem leads

* to a linear system, Ax=b, to be solved once in static problems and at each time step in

* time-dependent problems. There are two cases to be distinguished. The first case is where

* BOTH A and b are assembled in an FE manner, b possibly being composed of a volume integral

* plus a surface integral. The other case is where this is not true, for example where b = Mc+d,

* where the vector d and matrix M are assembled in an FE manner, but not the vector c.

*

* The Chaste PDE classes include ASSEMBLER classes, for setting up anything assembled in an

* FE manner, and SOLVER classes, for setting up linear systems. In the general case, solvers

* need to own assemblers for setting up each part of the linear system. However for the first

* case described above (in which both A and b in Ax=b are assembled), we can use the design

* where the solver IS AN assembler. We illustrate how to do this in the first tutorial.

*

* ### Writing solvers

*

* Let us write a solver for the coupled 2-unknown problem

*

* $$

* \begin{align*}

* \nabla^2 u + v &= f(x,y),\\\\\\

* \nabla^2 u + u &= g(x,y),

* \end{align*}

* $$

*

* where $f$ and $g$ are chosen so that,

* with zero-dirichlet boundary conditions, the solution is given by $u = \sin(\pi x)\sin(\pi x)$,

* $v = \sin(2\pi x)\sin(2\pi x)$.

*

* (As a brief aside, note that the solver we write will in fact work with general Dirichlet-Neumann

* boundary conditions, though the test will only provide all-Dirichlet boundary conditions. We

* save a discussion on general Dirichlet-Neumann boundary conditions for the second example.)

*

* Using linear basis functions, and a mesh with N nodes, the linear system that needs to be set up is

* of size 2N by 2N, and in block form is:

*

* ```

* [ K -M ] [U] = [b1]

* [ -M K ] [V] [b2]

* ```

*

* where `K` is the stiffness matrix, `M` the mass matrix, `U` the vector of nodal values

* of u, `V` the vector of nodal values of v, `b1` the vector with entries $\int f\phi_i dV$ ($i=1,\dots,N$)

* and `b2` has entries $\int g\phi_i dV$ (here $\phi_i$ are the linear basis functions).

*

* This is the linear system which we now write a solver to set up.

* {{< callout context="note" title="Note" icon="info-circle" >}}

* The main Chaste solvers assume a **STRIPED** data format, i.e. that the unknown vector

* is:

*

* `[U_1 V_1 U_2 V_2 .. U_n V_n]`, not `[U_1 U_2 .. U_n V_1 V_2 .. V_n]`.

*

* *We write down equations in block form as it makes things clearer, but have to remember that the code

* deals with STRIPED data structures.* Striping is used in the code for parallelisation reasons.

* {{< /callout >}}

*

* These are some basic includes as in the solving-PDEs tutorials

*/

#include <cxxtest/TestSuite.h>

#include "TetrahedralMesh.hpp"

#include "TrianglesMeshReader.hpp"

#include "BoundaryConditionsContainer.hpp"

#include "ConstBoundaryCondition.hpp"

#include "PetscSetupAndFinalize.hpp"

#include "TrianglesMeshWriter.hpp"

/* We need to include the following two classes if we are going to use a combination of

* (element_dim, space_dim, problem_dim) that isn't explicitly instantiated in

* `BoundaryConditionsContainer.cpp` (see the bottom of that file) (without these includes this test will

* fail to link).

*/

#include "BoundaryConditionsContainerImplementation.hpp"

#include "AbstractBoundaryConditionsContainerImplementation.hpp"

/* These two classes will be used in writing the solver */

#include "AbstractAssemblerSolverHybrid.hpp"

#include "AbstractStaticLinearPdeSolver.hpp"

/* We will solve a second problem, below, which will be time-dependent and will

* require the following class */

#include "AbstractDynamicLinearPdeSolver.hpp"

/* The linear system is Ax=b where A and b are FE assembled, so we can use the solver-is-an-assembler

* design. To construct our solver, we inherit from `AbstractAssemblerSolverHybrid` which links

* the solver clases to the assembler classes, and `AbstractStaticLinearPdeSolver`. (For time-dependent

* problems, the second parent would be `AbstractDynamicLinearPdeSolver`). Note the template

* parameter `PROBLEM_DIM` below, in this case it is 2 as there are two unknowns.

*/

class MyTwoVariablePdeSolver

};

/*

* That is the solver written. The usage is the same as see the PDE solvers described in the

* previous tutorials - have a look at the first test below.

*

* ### A solver of 3 parabolic equations

*

* Let us also write a solver for the following problem, which is composed of 3 parabolic PDEs

*

* $$

* \begin{align*}

* u_t &= \nabla^2 u + v,\\\\\\

* v_t &= \nabla^2 v + u + 2w,\\\\\\

* w_t &= \nabla^2 w + g(t,x,y),

* \end{align*}

* $$

*

* where $g(t,x,y) = t$ if $x>0.5$ and $0$ otherwise. This time we assume general

* Dirichlet-Neumann boundary conditions will be specified.

*

* The `AbstractAssemblerSolverHybrid` deals with the Dirichlet and Neumann boundary parts of the implementation,

* so, we, the writer of the solver, don't have to worry about this. The user has to realise though that they are

* specifying NATURAL Neumann BCs, which are whatever appears naturally in the weak form of the problem. In this case, natural

* Neumann BCs are specification of: $\partial u/\partial n = s_1, \partial v/\partial n = s_2, \partial w/ \partial n = s_3$,

* which coincide with usual Neumann BCs. However, suppose the last equation was $w_t = \nabla^2 w + \nabla . (D \nabla u)$,

* then the natural BCs would be:

* $\partial u/\partial n = s_1, \partial v/\partial n = s_2, \partial w/\partial n + (D\nabla u).n = s_3$.

*

* We need to choose a time-discretisation. Let us choose an implicit discretisation, ie

*

* ```

* (u^{n+1} - u^{n})/dt = Laplacian(u^{n+1}) + v^{n+1}

* (v^{n+1} - v^{n})/dt = Laplacian(v^{n+1}) + u^{n+1} + 2w^{n+1}

* (w^{n+1} - w^{n})/dt = Laplacian(w^{n+1}) + g(t^{n+1},x)

* ```

*

* Using linear basis functions, and a mesh with N nodes, the linear system that needs to be set up is

* of size 3N by 3N, and in block form is:

*

* ```

* [ M/dt+K -M 0 ] [U^{n+1}] = [b1] + [c1]

* [ -M M/dt+K -2M ] [V^{n+1}] [b2] + [c2]

* [ 0 0 M/dt+K ] [W^{n+1}] [b3] + [c3]

* ```

*

* where `K` is the stiffness matrix, `M` the mass matrix, `U^n` the vector of nodal values

* of u at time t_n, etc, `b1` has entries `integral( (u^n/dt)phi_i dV )`, and similarly for

* `b2` and `b3`. Writing the Neumann boundary conditions for

* u as `du/dn = s1(x,y)` on Gamma, a subset of the boundary, then `c1` has entries

* `integral_over_Gamma (s1*phi_i dS)`, and similarly for `c2` and `c3`.

*

* Let us create a solver for this linear system, which will be written in a way in which the RHS

* vector is assembled in an FE manner, so that the solver-is-an-assembler design can be used.

* Note that this solver inherits from `AbstractDynamicLinearPdeSolver` and PROBLEM_DIM is now equal

* to 3. We don't have to worry about setting up [c1 c2 c3] (we just need to take in a `BoundaryConditionsContainer`

* and the parent will use it in assembling this vector). We do however have to tell it

* how to assemble the volume integral part of the RHS vector, and the LHS matrix.

*/

class ThreeParabolicPdesSolver

};

/* Now the tests using the two solvers */

class TestWritingPdeSolversTutorial : public CxxTest::TestSuite

};

/*

*

* **Visualisation:** To visualise in matlab/octave, you can load the node file,

* and then the data files. However, the node file needs to be edited to remove any

* comment lines (lines starting with `#`) and the header line (which says how

* many nodes there are). (The little matlab script [`anim/matlab/read_chaste_node_file.m`](https://github.com/Chaste/Chaste/blob/develop/anim/matlab/read_chaste_node_file.m)

* can be used to do this, though it is not claimed to be very robust). As an

* example matlab visualisation script, the following, if run from the output

* directory, plots w:

*

* ```matlab

* pos = read_chaste_node_file('mesh.node');

* for i=0:200

* file = ['txt_output/results_Variable_2_',num2str(i),'.txt'];

* w = load(file);

* plot3(pos(:,1),pos(:,2),w,'*');

* pause;

* end;

* ```

*

*/

#endif // TESTWRITINGPDESOLVERSTUTORIAL_HPP_