PyEFVLib 1.0.12

A library that facilitates the implementation of the Element-based Finite Volumes Method (EbFVM) to solve partial differential equations (PDEs)

PyEFVLib

PyEFVLib is a Python library that provides tools to support the solution of systems of partial differential equations (PDEs) using the Element-based Finite Volumes Method (EbFVM)

Summary

1. Installation
2. Numerical Method
3. PyEFVLib Classes - Geometrical Entities
4. PyEFVLib Classes - Non-Geometrical Entities
5. Tutorial
6. 3rd-party Softwares
7. 2D vs 3D considerations
8. bellbird

1. Installation

Install PyEFVLib using pip, by typing in the command line:

pip install PyEFVLib

2. Numerical Method

2.1 - Partial Differential Equations

• A differential equation is an equation that relates one or more functions and their derivatives.
• A partial differential equation (PDE) is an equation that relates one or more functions dependent on two or more variables, to the partial derivatives of these functions in relation to the independent variables.
• In physics, differential equations relate physical properties (e.g., pressure, temperature, displacement, turbulent kinetic energy, speed) to their spatial and temporal derivatives, and model natural and everyday phenomena.

---

2.2 - Example: Integrating a PDE in space

Take the heat equation as an example

Integrating in the control volume (CV)

And applying the divergence theorem

---

2.3 - Domain Discretization

Previously it was shown that a differential equation in the differential form can be integrated in any control volume, but to achieve a detailed solution along the domain of interest, the finite volume method discretizes the domain in several control volumes forming a mesh. The more refined the mesh, the better the solution to the problem will be.

Note: The black bordered triangles shown in the figure are actually the elements, the control volumes are outlined by the blue lines.

---

2.4 - Field Approximations

As it is the nature of PDEs, the equations involve spatial derivatives, but since we are storing a finite number of variables, we need a way to approximate the field continuously to evaluate the spatial derivatives, and the EbFVM does this using the shape functions.

We'll start by evaluating the field within an element of our mesh. It'll be expressed by a linear combination of the property values at the vertices of the element, and the weight of the property at each vertex is given by the shape functions at each point inside the element, as in the figure below:

And as we know the entire field inside the element, we can evaluate its partial derivatives, as shown below:

---

2.5 - Problem Solution

Previously, the differential equations expressed the values of the unknowns as a function of their partial derivatives, but as we can express the partial derivatives as a function of the values of the field at the neighboring vertices, so the solution of the problem becomes a system of algebraic equations.

And in the case of a linear system of equations

---

3. PyEFVLib Classes - Geometrical Entities

3.1 - Grid

Attribute	Description
vertices	a list of the grid's vertices
elements	a list of the grid's elements
regions	a list of the grid's regions
boundaries	a list of the grid's boundaries
dimension	the dimension of the mesh (2 or 3)
gridData	a data structure that holds the information provided by mesh file

---

3.2 - Vertex

Attribute	Description
x	vertex x coordinate
y	vertex y coordinate
z	vertex z coordinate
handle	the vertex index in the grid
elements	a list of all elements to which vertex belongs
volume	the volume of the control volume at which vertex is centered
coordinates	a numpy array with the vertex coordinates

Method	Parameters	Return type	Description
getLocal	element	int	returns the local index of the vertex within the element
getInnerFaces		list[InnerFace]	returns the innerFaces next to vertex
getSubElementVolume	element	float	returns the volume of the section of element that is inside vertex control volume

Vertex is a child class of the class Point

It is possible to do operations like vertex1 + vertex2, a * vertex, ...

---

3.3 - Region

Region is a class used to perform simulations where there are multiple regions with different properties, so the properties are specific to each region. If the domain is homogeneous (has constant properties) one single region will be created containing all the elements. The region names are set in the mesh file.

Attribute	Description
name	the name of the region
handle	the region index in the grid
elements	a list of all elements belonging to region

---

3.4 - Element

Attribute	Description
handle	the element index in the grid
region	the region to which element belongs
vertices	a list of the vertices of element
innerFaces	a list of the element's innerFaces
outerFaces	a list of the element's outerFaces
faces	a list of the element's innerFaces and outerFaces
shape	the element Shape
volume	the sum of all subElementVolumes (control volume section inside element)
subElementVolumes	the list with each vertex subElementVolume
centroid	the element centroid
numberOfVertices	the number of vertices of element

Method	Parameters	Return type	Description
getTransposedJacobian	shapeFunctionDerivatives	numpy.array[d,m]	returns the transposed of the jacobian of shapeFunctionDerivatives
getGlobalDerivatives	derivatives	numpy.array[d,d]	returns the derivatives of the vertices shape functions with respect to x,y,z

Where d is the grid dimension and m is the number of vertices in the element

---

3.5 - InnerFace

The innerFaces are inside the element, but they surround the control volumes. The reason we need the shape functions and its derivatives evaluated at the inner faces centroids is because whenever there's a surface integral in our differential equation (over the control surface), it'll be summed over "integration points" over the control surface, which happen to be at the centroids of the inner faces.

Attribute	Description
handle	the innerFace index in the grid
local	the innerFace local index in the element
element	the innerFace element
area	the innerFace area (Point object)
centroid	the innerFace centroid (Point object)
globalDerivatives	a numpy.array[d,m] containing the gradient operator. Multiply by an m-dimensional array containing the property at the vertices of its element to get the gradient of the property at the innerFace

Method	Parameters	Return type	Description
getShapeFunctions	-	numpy.array[d,m]	The shape functions of the element's vertices' shapeFunctions evaluated at the innerFace
getNeighborVertices	-	(Vertex,Vertex)	The surrounding vertices of innerFace
getNeighborVerticesLocals	-	(int, int)	The surrounding vertices' local indices in the element of innerFace
getNeighborVerticesHandles	-	(int, int)	The surrounding vertices' handles of innerFace

---

3.6 - Boundary

Attribute	Description
name	it used to apply the boundary conditions on it
handle	the boundary index on the grid
vertices	a list with the vertices that make the boundary
facets	a list with the facets that make the boundary

---

3.7 - Facet

Attribute	Description
handle	the facet index on the grid
boundary	the boundary to which facet belongs
element	the element to which facet belongs
vertices	the vertices that lie on facet
outerFaces	the outerFaces that make facet
area	facet's area (Point object)
centroid	facet's centroid (Point object)
elementLocalIndex	its local index on element
boundaryLocalIndex	its local index on boundary

---

3.8 - OuterFace

Attribute	Description
local	local index on facet
vertex	neighboring vertex
facet	facet to which outerFace belongs
centroid	its centroid (Point object)
area	its area (Point object)

Method	Parameters	Return type	Description
getShapeFunctions	-	numpy.array[d,m]	The shape functions of the element's vertices' shapeFunctions evaluated at the outerFace

---

3.9 - Shape

The Shape subclasses are:

• Triangle
• Quadrilateral
• Tetrahedron
• Hexahedron
• Prism
• Pyramid

All of them share the same attributes and methods:

Attribute	Description
dimension	the shape dimension
numberOfInnerFaces	the number of inner faces
numberOfFacets	the number of facets
subelementTransformedVolumes	the volumes of the subelements in the transformed space
innerFaceShapeFunctionValues	the shape function values evaluated at the innerFace centroid
innerFaceShapeFunctionDerivatives	the shape function derivatives evaluated at the innerFace centroid
innerFaceNeighborVertices	the element local indices of the vertices that surround each innerFace
subelementShapeFunctionValues	the shape function values evaluated at the subelement centroid
subelementShapeFunctionDerivatives	the shape function derivatives evaluated at the subelement centroid
facetVerticesIndexes	the element local indices of the group of vertices that make each facet
outerFaceShapeFunctionValues	the shape function values evaluated at the outerFace centroid
vertexShapeFunctionDerivatives	the shape function derivatives evaluated at the vertices

---

3.10 - BoundaryConditions

The two BoundaryConditions classes are

• DirichletBoundaryCondition
• NeumannBoundaryCondition

Attribute	Description
handle	the boundaryCondition index on the grid
value	the prescribed condition value
expression	True if value is a string and the expression must be parsed, otherwise False
boundary	the boundary to which the condition is being applied

Method	Parameters	Return type	Description
getValue	index, time=0.0	float	if expression==False, returns value. Otherwise, will parse the value string

If expression==True, the value string is an expression that will be parsed using the built-in method eval. It may contain the variables x, y and z, which are equal respectively to grid.vertices[index].x, grid.vertices[index].y and grid.vertices[index].z, as well as t, which is equal to the argument time. The following functions/constants can also be used: pi, sin, cos, tan, arcsin, arccos, arctan, sinh, cosh, tanh, arcsinh, arccosh, arctanh, sqrt, e , log, exp, inf, mod, floor (which have been imported from the module numpy).

The Dirichlet boundary condition means that the unknown variable is being prescribed at the boundary, and the Neumann boundary condition means that the unknown variable directional gradient (in the direction of the normal vector to the boundary) is being prescribed.

---

4. PyEFVLib Classes - Non-Geometrical Entities

Besides grid entities generation that facilitates the solution of PDEs, PyEFVLib also offers some I/O classes, some of which are integrated with other I/O libraries such as meshio.

4.1 - ProblemData

The ProblemData class purpose is to provide the solver every information needed to solve the problem, in such a way that one can write a code in the following way:

def solver(problemData):
  ...
  
if __name__ == "__main__":
  problemData = PyEFVLib.ProblemData(...)
  solution = solver(problemData)

This way in order to keep everything simple ProblemData only needs 5 arguments to initialize

Arguments	Arg Type	Description
meshFilePath	string	provides the path to the input mesh
outputFilePath	string	provides the directory where the results will be written
numericalSettings	NumericalSettings	provides all the numerical method details
propertyData	PropertyData	provides all the property info of the domain
boundaryConditions	BoundaryConditions	provides all the boundary conditions of the problem

• 4.1.1 - NumericalSettings

Arguments	Arg Type	Default	Description
timeStep	float	-	step in time to be increased
finalTime	float	inf	when the simulation reaches this time it'll stop
tolerance	float	1e-4	can be a steady state criterion or an iterative tolerance (it is up to who writes the solver and can be used in both ways)
maxNumberOfIterations	int	1000	when the simulation reaches this number of iterations it'll stop

• 4.1.2 - PropertyData

Usage:

numericalSettings = PyEFVLib.NumericalSettings({
  "RegionName1": {
    "property1": value1,
    "property2": value2,
    "property3": value3,
  },
  "RegionName2": {
    "property1": value4,
    "property2": value5,
    "property3": value6,
  },
})

• 4.1.3 - BoundaryConditions

Usage:

boundaryConditions = PyEFVLib.BoundaryConditions({
  "variable1": {
    "InitialValue": initialValue1,
    "boundaryName1": { "condition": PyEFVLib.Dirichlet, "type": PyEFVLib.Constant, "value": 0.0 },
    "boundaryName2": { "condition": PyEFVLib.Dirichlet, "type": PyEFVLib.Variable, "value": "100*sin(x*t)" },
    "boundaryName3": { "condition": PyEFVLib.Neumann, "type": PyEFVLib.Constant, "value": 100.0 },
  },
  "variable2": {
    "InitialValue": initialValue2,
    "boundaryName1": { "condition": PyEFVLib.Neumann, "type": PyEFVLib.Variable, "value": "0.0 if y<0.5 else 50.0" },
    "boundaryName2": { "condition": PyEFVLib.Dirichlet, "type": PyEFVLib.Constant, "value": -50.0 },
    "boundaryName3": { "condition": PyEFVLib.Neumann, "type": PyEFVLib.Constant, "value": 0.0 },
  },
})

---

4.2 - Readers

The current readers are:

• MSHReader
• XDMFReader

(Soon - MeshioReader)

The only argument that the readers take is the mesh path. To initialize a grid object from a reader do as such:

grid = PyEFVLib.Grid( MSHReader(filePath).getData() )

or

grid = PyEFVLib.read(filePath)

To create your own reader, all it is needed is to end your reader by creating a GridData object (which only takes the mesh file path as argument), and populate it by using the methods: setDimension, setVertices, setElementConnectivity, setRegions, setBoundaries, setShapes

---

4.3 - Savers

The current savers are:

• CsvSaver (saves the results in a .csv spreadsheet-like file)
• VtuSaver (saves the results in only one moment in a .vtu file)
• VtmSaver (saves the results in multiple .vtm files from all time steps)
• MeshioSaver (mostly saves the results in one moment, except by xdmf) [recommended]

MeshioSaver supports the following extensions: msh, mdpa, ply, stl, vtk, vtu, xdmf, xmf, h5m, med, inp, mesh, meshb, bdf, fem, nas, obj, off, post, post.gz, dato, dato.gz, su2, svg, dat, tec, ugrid, wkt. Works best with xdmf

Arguments	Arg Type	Description
grid	Grid	the grid object
outputFilePath	str	provided by problemData.outputFilePath
libraryPath	str	provided by problemData.libraryPath
extension	str	exclusive of MeshioSaver, insert one of the listed above

Method	Parameters	Description
save	variableName, field, currentTime	saves the field
finalize	-	MUST be called at the end of the simulation in order for the results to be saved

---

5. Tutorial

In order to demonstrate the library classes we'll solve the heat transfer equation. Let's remember the equations:

We'll start by importing the library and initializing some variables

import PyEFVLib
import numpy as np

def heatTransfer(problemData):
    propertyData     = problemData.propertyData
    timeStep         = problemData.timeStep
    grid             = problemData.grid
    numberOfVertices = problemData.grid.numberOfVertices
    dimension        = problemData.grid.dimension
  
    saver = PyEFVLib.MeshioSaver(grid, problemData.outputFilePath, problemData.libraryPath, extension="xdmf")
  
    temperatureField = np.repeat(0.0, numberOfVertices)
    oldTemperatureField = problemData.initialValues["temperature"].copy()

The first 5 lines where we call some problemData variables are optional, but since they are used very often, it is handy to declare them like that. The saver must be initialized before starting to save the results, so here is a good place to do it. The temperatureField is our unknown to be found, and the oldTemperatureField will be updated after each increment in time, and will be used to evaluate the temporal derivative.

Notice that in our discretized equation the left-hand side contains only terms involving the temperature (which are going to be added to the matrix) and the right-hand side contains only terms independent of the temperature(and are going to be added to the independent vector). The linear system matrix will have N rows and N columns, where the i-th row corresponds to the equation applied to the i-th control volume (or vertex), and the columns correspond to the temperature in the i-th vertex. So, to implement the matrix of our linear system we'll do the following:

    ...
    oldTemperatureField = problemData.initialValues["temperature"].copy()
    
    def assembleMatrix():
        matrix = np.zeros((numberOfVertices, numberOfVertices))
    
        for region in grid.regions:
            k     = propertyData.get(region.handle, "conductivity")
            q     = propertyData.get(region.handle, "heatGeneration")
            rho = propertyData.get(region.handle, "density")
            cp    = propertyData.get(region.handle, "specificHeat")
            
            for element in region.elements:
                # Accumulation term (rho * cp * T/Δt)
                for vertex in element.vertices:
                    matrix[vertex.handle][vertex.handle] += vertex.getSubElementVolume(element) * rho * cp / timeStep
                    
                # Diffusion term ( -k*grad(T) Δs <-> -k * Δs^T * B_ip * T_elem )
                for innerFace in element.innerFaces:
                    area = innerFace.area.getCoordinates()[:dimension]
                    backwardsVertexHandle, forwardVertexHandle = innerFace.getNeighborVerticesHandles()
                    
                    coefficients = -k * np.matmul(area.T, innerFace.globalDerivatives)
                    for local, vertex in enumerate(element.vertices):
                        matrix[backwardsVertexHandle][vertex.handle] += coefficients[local]
                        matrix[forwardVertexHandle][vertex.handle] -= coefficients[local]

Remember that the region object contains the elements with the same properties, that's why we start the assembleMatrix function looping over the regions and calling the properties at that region. Then we loop over the elements, since the innerFaces belong inside the elements, and for the vertices loop we need the subElementVolume.

The term comes from the i-th equation relating the i-th vertex temperature, that's why we add the term to matrix[vertex.handle][vertex.handle]

Also note that from our surface integral over the control surface, we got the summation over the centroids of the innerFaces (or integration points, denoted in the equations as ip). So the loop over the innerFaces of the element will add this term to our linear system of equations. The (or innerFace.area) is the area vector of the innerFaces, which modulus is equal to the area of the innerFace and points in the outwards normal direction to the innerFace. So in our summation we add coefficient to backwardsVertex because innerFace.area is pointing outwards from its control volume, and we subtract it from forwardVertex because -innerFace.area is pointing outwards from its control volume.

Since innerFace.area is a Point object (and not a numpy.array) we convert it into a 3x1 numpy.array using the getCoordinates method, and because innerFace.globalDerivatives has dimension rows and element.numberOfVertices columns, we slice it into a dimensionx1 numpy.array (removing the z coordinate if the simulation is 2D).

The term is called gradient operator, and it is stored in innerFace.globalDerivatives. It's defined as

and it has the following property: .

Now we need to turn our attention to the boundary conditions. Previously we simplified our equation by discretizing it only in the interior of the domain, without making considerations about its boundaries. Note that we can express Control Surface = Control Surface (inside the domain) + Control Surface (where the Neumann boundary condition is applied) + Control Surface (where the Dirichlet boundary condition is applied).

So in the step we could have made the following substitution:

Where

And

The Dirichlet boundary condition is very simple, if it is applied to the i-th vertex, all the elements of the i-th row of the matrix must be 0, except for the i-th element of the i-th row, which is equal to 1. Meanwhile, the i-th element of the independent vector must be the prescribed temperature.

That being said, the way we implement the boundary conditions in our matrix is the following

(still in the assembleMatrix function)

...
def heatTransfer(problemData):
    ...
    def assembleMatrix():
        ...
                        matrix[forwardVertexHandle][vertex.handle] += k * np.matmul(area.T, innerFace.globalDerivatives)
                        
        # Dirichlet Boundary Condition
        for boundaryCondition in problemData.dirichletBoundaries["temperature"]:
            for vertex in boundaryCondition.boundary.vertices:
                matrix[vertex.handle] = np.zeros(numberOfVertices)
                matrix[vertex.handle][vertex.handle] = 1.0
        
        return matrix

Now we need an assembleIndependentVector function:

...
def heatTransfer(problemData):
    ...
    def assembleIndependentVector():
        independent = np.zeros(numberOfVertices)
    
        for region in grid.regions:
            k   = propertyData.get(region.handle, "conductivity")
            q   = propertyData.get(region.handle, "heatGeneration")
            rho = propertyData.get(region.handle, "density")
            cp  = propertyData.get(region.handle, "specificHeat")
            
            for element in region.elements:
                for vertex in element.vertices:
                    # Accumulation term (rho * cp * T_old/Δt)
                    independent[vertex.handle] += vertex.getSubElementVolume(element) * rho * cp * oldTemperatureField[vertex.handle] / timeStep
                    
                    # Generation term (q''')
                    independent[vertex.handle] += vertex.getSubElementVolume(element) * q
                    
            # Neumann Boundary Condition
            for boundaryCondition in problemData.neumannBoundaries["temperature"]:
                for facet in boundaryCondition.boundary.facets:
                    for outerFace in facet.outerFaces:
                        independent[outerFace.vertex.handle] += boundaryCondition.getValue(outerFace.vertex.handle, currentTime) * np.linalg.norm(outerFace.area.getCoordinates())
                        
            # Dirichlet Boundary Condition
            for boundaryCondition in problemData.dirichletBoundaries["temperature"]:
                for vertex in boundaryCondition.boundary.vertices:
                    independent[vertex.handle] = boundaryCondition.getValue(vertex.handle, currentTime)
                    
        return independent

Because the assembleMatrix function and the boundary condition details have already been explained, this function is more straightforward.

The accumulation term is added to the independent vector, and it calls the oldTemperatureField vector, which must be updated after each step in time.

Both the Neumann boundary condition and the Dirichlet boundary condition are applied to the facets of the boundaries, which are in a way different from the vertices. However, since the Dirichlet boundary condition prescribes directly the temperature (or whichever property) at the vertices of the boundary, the temperature will be set at those points, that's why it is applied after the Neumann boundary condition, because otherwise the temperature would be mistakenly set.

On the boundaryCondition.getValue(outerFace.vertex.handle), remember that if the boundaryCondition is of the type PyEFVLib.Constant, it will return the prescribed boundaryCondition anyways, however if it is of the type PyEFVLib.Variable, it can depend on the x,y and z coordinates of the vertex and also of the currentTime. So we must pass the handle of the outerFace.vertex so it can check the x,y,z coordinates of it.

Next we start the loop to solve our problem:

    steadyStateCriterion = problemData.tolerance
    difference = 2*steadyStateCriterion
    stepInTime = 0
    currentTime = 0.0
    convergedToSteadyState = False
    
    matrix = assembleMatrix()
    while not convergedToSteadyState:
        independent = assembleIndependentVector()
        temperatureField = np.linalg.solve(matrix, independent)
        difference = max( abs(temperatureField - oldTemperatureField) )
        
        oldTemperatureField = temperatureField.copy()
        currentTime += timeStep
        stepInTime += 1
        
        saver.save("temperatureField", temperatureField, currentTime)
        print("{:>9}        {:>14.2e}     {:>14.2e}     {:>14.2e}".format(stepInTime, currentTime, timeStep, difference))
        
        convergedToSteadyState = ( difference <= steadyStateCriterion ) or ( currentTime >= problemData.finalTime ) or ( stepInTime >= problemData.maxNumberOfIterations )
            
    saver.finalize()

Here we declare some handy variables at the start of the loop, and notice that we assemble the matrix before starting the loop. That's because none of the terms added in matrix change over time. If we wanted for example to have a variable timeStep, then we'd have to assemble the matrix inside the transient loop every instant of time. In this example, other detail could have been changed for the simulation to perform better. Instead of solving the linear system of equations every single step in time, we could have inverted the matrix before the transient loop (right after assembling the matrix) and simply multiplied the inverted matrix by the independent vector. Here it has not been done for didactic purposes.

The difference variable is meant to track how the temperatureField stop changing after a while, and to stop the simulation after the changes are negligible. The .copy() method is needed to update oldTemperatureField because of how the mutable list class behaves. After solving the problem it is needed to call the saver.finalize() method, to write the output file with the results.

Finally, we need to call our heatTransfer function, like so:

if __name__ == "__main__":
    problemData = PyEFVLib.ProblemData(
        meshFilePath = "mesh.msh",
        outputFilePath = "results/",
        numericalSettings = PyEFVLib.NumericalSettings(timeStep = 1000.0, tolerance = 1e-4, maxNumberOfIterations = 500),
        propertyData = PyEFVLib.PropertyData({
            "Body": {
                "conductivity"   : 22.0,   # [W/m.K]
                "density"        : 8960.0, # [kg/m3]
                "specificHeat"   : 377.0,  # [J/kg.K]
                "heatGeneration" : 5000.0, # [W/m3]
            },
        }),
        boundaryConditions = PyEFVLib.BoundaryConditions({
            "temperature": {
                "InitialValue": 300.0,
                "West": { "condition" : PyEFVLib.Dirichlet, "type" : PyEFVLib.Constant, "value" : 300.0 },
                "East": { "condition" : PyEFVLib.Dirichlet, "type" : PyEFVLib.Constant, "value" : 400.0 },
                "South": { "condition" : PyEFVLib.Neumann, "type" : PyEFVLib.Constant, "value" : 0.0 },
                "North": { "condition" : PyEFVLib.Neumann, "type" : PyEFVLib.Constant, "value" : 0.0 },                
            }
        }),
    )
    
    heatTransfer(problemData)

---

6. 3rd-party Softwares

6.1 - Gmsh

Notice that in the workflow presented, the mesh is simply imported by passing its path, but it needs to be generated. In order to do so, the recommended software is Gmsh, which is an open-source mesh generator with built-in CAD tools. The only (and minor) downside of the mesh importing at the moment is that the mesh must be from the 2.2 version of .msh. But that's very simple. Just download the latest version of Gmsh, and follow the step-by-step:

1. Open the Help menu
2. Choose Current Options and Workspace
3. Search 'version' in the search bar, or scroll down until you find the parameter 'Mesh.MshFileVersion'
4. Double click it, type 2.2, and click Ok
5. Optionally open the File menu, and choose the option 'Save Options as Default', so you don't have to go through this step-by-step every time

6.2 - Paraview

After solving the simulation, you'll want to view and analyze the results you generated. One option is to save it as a .csv file (using the PyEFVLib.CsvSaver class) and open it using pandas, and you could visualize it using matplotlib. But if you want to view your results, a much better way of doing this is with Paraview.

Paraview is an open-source data analysis and visualization application, and by saving the results in formats such as .xdmf (recommended), .vtu, .vtm, you can visualize it using paraview.

---

7. 2D vs 3D considerations

One great advantage of using the EbFVM, is that the difference between 2D and 3D simulations is very little. The major changes are in the field and gradient approximations using different elements and therefore different shape functions. However, PyEFVLib handles it without the user having to pay attention. For example in our heat transfer application show as an example in section 5, the only change one need to make to simulate a 3D domain is simply to change the mesh file path (and generate a 3D mesh), and add the boundary conditions according to the new mesh.

---

8. bellbird

Although PyEFVLib makes the use of EbFVM quite simple for the user, it is still required for the user to know the numerical method, how to discretize the equations and the full structure of the library. And when solving bigger systems of differential equations the coupling of the equations in the linear system can become a quite tricky, and can demand more attention to little details, while also requiring a lot more lines of code. For that reason, the bellbird library makes it way more simple for the user to solve systems of PDEs. The user informs the equations as string, and the library parses, integrates and discretizes them, and writes a script very similar to one the shown in the tutorial using PyEFVLib.

Then, if the user wants, they can change the code, but with a huge headstart.

At the moment, the library is still under development, but it can already solve some systems of equations and offer a solid basis for anyone who wants to solve their own problem. Check out the usage for the problem we solved earlier:

import bellbird

model = bellbird.Model(
	name = "Heat Transfer",
	equationsStr = ["k * div( grad(T) ) + q = rho * cp * d/dt(T)"],
	variables   = ["T"], # temperature [K]
	properties = {
		"Body":{
			"k" : 22.0,	# [W/m.K]	- conductivity
			"rho" : 8960.0,	# [kg/m3]	- density
			"cp" : 377.0,	# [J/kg.K]	- specific heat
			"q" : 5000.0,	# [W/m3]	- heat generation
		},
	},
	boundaryConditions = [
		bellbird.InitialCondition("T", 0.0),
		bellbird.BoundaryCondition("T", bellbird.Dirichlet, "West", 300.0),
		bellbird.BoundaryCondition("T", bellbird.Dirichlet, "East", 400.0),
		bellbird.BoundaryCondition("T", bellbird.Neumann, "South", 0.0),
		bellbird.BoundaryCondition("T", bellbird.Neumann, "North", 0.0),
	],
	meshPath = "mesh.msh",
	timeStep = 1000.0,
	tolerance = 1e-4,
	maxNumberOfIterations = 500,
	sparse = False,
)

model.run()