# Physical Solvers

Amongst the main components required in a topology optimization dedicated plateform, there is the so-called physical solver brick. Its purpose it to provide a numerical approximation of the solution to a partial differential equation corresponding to the problem at hand, given a set of boundary condition and on a given mesh, ergo a discretisation of a domain \(\Omega\).

The `OpenPisco.PhysicalSolvers.SolverBase`

class aims at providing a base structure and a common interface to the physical solvers usable within the platform.

## Auxiliary field/scalar

When solving a physical problem, at the very least, we expect the physical solver to compute the solution of the PDE at hand (be it a displacement, a temperature or something else). However, depending on the use case, we may want to rely on the physical solver to compute in the same time other physical quantities. This is the issue addressed by the auxiliary concept.

Note that the current solvers available in OpenPisco do not have the exact same capabilities, ergo some may be able to compute natively on their own (with or without a little help) a given physical quantity while the others may not.

The SolverBase class contains a dictionary named `OpenPisco.PhysicalSolvers.SolverBase.SolverBase.auxiliaryFieldGeneration`

to store whether a solver is able to provide if required a given field evaluated on different support (namely on nodes, centroids or integration points for instance). In a derived class, such an attribute is redefined according to the solver capabilities.

Thus

To check whether a solver is able to compute a named quantity on a given support, one can simply use

`OpenPisco.PhysicalSolvers.SolverBase.SolverBase.CanSupplyAuxiliaryField()`

Before solving the problem, one can ask the problem to generate a named quantity on a given support through

`OpenPisco.PhysicalSolvers.SolverBase.SolverBase.SetAuxiliaryFieldGeneration()`

. These quantities are to be computed when the problem is actually solved only.After Solving the physical problem, one can retrieve any named quantity on a given support after the problem has been solved but that were asked prior to the computation through

`OpenPisco.PhysicalSolvers.SolverBase.SolverBase.GetAuxiliaryField()`

A similar logic is used for scalar quantities.

## FieldsNames

All the fields and scalars mentionned are defined in `OpenPisco.PhysicalSolvers.FieldsNames`

. The idea is to establish a
mapping between each field/scalar variable and the actual field/scalar name (a string) used to denote a given field in physical solvers.

For instance, this is the name expected when one tries to retrieve a field/scalar from a file generated by a solver through the solver interface.

## SolveByLevelSet

In the context of level set-based topology optimization, two settings are adressed:

The conform case, or body-fitted approach, where the domain \(\Omega\) is simply bounded by the 0-th isusurface of the level set function. It can be extracted with

`OpenPisco.PhysicalSolvers.SolverBase.SolverBase.ExtractInteriorMesh()`

The fix mesh setting, or non conform, where we do not have an explicit mesh of the domain \(\Omega\) and therefore have to rely on the Ersatz material approximation (the default value for this pseudo-material stiffness is given in

`OpenPisco.PhysicalSolvers.SolverBase.SolverBase.eVoid`

)

Depending of the type of level set-based topology optimization performed, the mesh given to the physical solver (so-called `OpenPisco.PhysicalSolvers.SolverBase.SolverBase.cleanMesh`

) must be defined accordingly in the derived classes. Note that it does not mean all physical analysis supported by the physical solvers available are supposed to handle both cases.
The last step is to call the actual physical solver to solve the physical problem.

As of now, 4 interfaces with finite element solvers are available (three of them relies on external tools while the last is a native solver). Even though they do share a common interface, the core principles behind their respective implementation is different. The Macro & template files related to three of them are in the ExternalTools module.

## Physical solvers available

### External tools

Albeit slightly different regarding their core inner working, these solvers still share a common interface.

#### Code_Aster

Code_Aster 1 is an open-source finite element solver developed by EDF, for the analysis of structures and thermomechanics. Its core implementation is in Fortran and is also endowed as a domain specific language as a commande language. It can handle many behavior laws (elasticity, thermoelasticity, plasticity, fracture, hyperelasticity), finite elements (continuous 2D/3D, structural elements, kinematic couplings) and types of analysis (linear/nonlinear static, quasi-static, modal, dynamic on a modal basis or direct). It is also possible to natively compute quantities at Gauss points or at nodes.

The interfacing with OpenPisco is done with the following entities:

The class

`OpenPisco.ExternalTools.Aster.AsterInterface.AsterInterface`

, whose main responsabilities is toSelect and retrieve the proper study among the template available based on the type of physical analysis considered, using

`OpenPisco.ExternalTools.Aster.AsterInterface.AsterInterface.GetStudyFilesFromPath()`

. The studies are all located in`OpenPisco.ExternalTools.Aster`

whereThe *.comm files are scripts using Code_Aster domain specific high-level language to define and solve a physical problem

The *.export are the actual file given to the executable. It contains input/output filenames required for/produce by the physical analysis, including the *.comm files ordered adequately for the analysis to run properly

Call the executable to solve the physical problem at hand with

`OpenPisco.ExternalTools.Aster.AsterInterface.AsterInterface.ExecuteAster()`

Read the files generated by Code_Aster in order, for instance, to retrieve fields through

`OpenPisco.ExternalTools.Aster.AsterInterface.AsterInterface.GetField()`

Classes derived from

`OpenPisco.ExternalTools.Aster.AsterSolverBase.AsterSolverBase`

, for instance in`OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver`

, where each of them are associated to a specific physical analysis.

Their role is to:

Generate a proper python script to be used by the studies (*.comm) with

`OpenPisco.ExternalTools.Aster.AsterSolverBase.AsterSolverBase.WriteParametersInput()`

. Amongst others, and depending on the type of analysis, it can be employed to provide information to the solver regardingmaterial properties

boundary regions

boundary conditions

Specify the auxiliary fields/scalars to be generated with

`OpenPisco.ExternalTools.Aster.AsterSolverBase.AsterSolverBase.WriteAuxiliaryInputs()`

Solve the physical problem with

`OpenPisco.ExternalTools.Aster.AsterSolverBase.AsterSolverBase.Solve()`

#### Freefem++

Freefem++ 2 is a high level finite element software using a C++ idiom domain specific language. It aims to ease the setting of a physical problem by relying on its variational formulation through this idiom.

The interfacing with OpenPisco consists in:

The class

`OpenPisco.ExternalTools.FreeFem.FreeFemInterface.FreeFemInterface`

, whose main responsabilities is toCall the executable to solve the physical problem at hand with

`OpenPisco.ExternalTools.FreeFem.FreeFemInterface.FreeFemInterface.ExecuteFreefem()`

Read the files generated by FreeFem in order, for instance, to retrieve fields through

`OpenPisco.ExternalTools.FreeFem.FreeFemInterface.FreeFemInterface.GetNodalField()`

The class

`OpenPisco.ExternalTools.FreeFem.FreeFemProblemWriter.FreeFemProblemWriter`

, whose main responsabilities is to generate the Freefem scripts by combining the physical problem information in`OpenPisco.ExternalTools.FreeFem.FreeFemProblemWriter.FreeFemProblemWriter.WriteFreefemProblem()`

, such asmaterial properties

boundary regions

boundary conditions

#### Zset

Zset 3 is…To be described

### Internal tools

Muscat solver 4 is the internal finite element solver available within OpenPisco. It is implemented in Python, and Cython/C++ for performance critical routines and relies on the problem variational formulation. The interfacing with OpenPisco consist in a single module `OpenPisco.PhysicalSolvers.UnstructuredFEAGenericProblem`

.

Note that, unlike the external tools:

There is no need to performed any output operation like writing the mesh for instance, since OpenPisco already use Muscat routine to handle the mesh

There is no executable involved

## Physical analysis available

We provide now a list of the physical analysis supported by the platform, it is mean to be exhaustive. Let \(\Omega \subset \mathbb{R}^3\) such that \(\partial \Omega = \Gamma_D \cup \Gamma_N \cup \Gamma\).

### Linear elasticity analysis

We denote by \(u\) the solution in displacement of the following problem

where:

\(\sigma\) is the stress tensor satisfying the Hooke’s law

\(e(u)\) is the strain tensor in the small deformation framework

\(\lambda\) and \(\mu\) are the Lamé coefficient and can also be expressed with respect to \(E\), the Young modulus, and \(\kappa\), the Poisson ratio,

\(g\) is the density of surfaces forces acting on \(\Gamma_N\)

\(f\) is the density of volumic forces acting on \(\Omega\)

the body is clamped on \(\Gamma_D\)

Associated implementation:

`OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterStatic`

`OpenPisco.PhysicalSolvers.GeneralZSetSolverMeca.GeneralZSetPhysicalProblemMeca`

`OpenPisco.PhysicalSolvers.StaticFreeFemSolverMeca.StaticFreeFemSolverMeca`

`OpenPisco.PhysicalSolvers.UnstructuredFEAGenericSolver.UnstructuredFEAGenericProblem`

### Thermal analysis

We are looking for the temperature \(T\), solution of the following problem

where:

\(S\) is the thermal source

\(T_{ext}\) is the emperature of the fluid in contact with the structure

\(T_d\) is the prescribed temperature

\(h\) is the convection coefficient

\(\rho\) is the material density

\(C_p\) is the thermal capacity

\(\lambda_c\) is the thermal conductivity

\(T_{init}\) is the initial temperature

the temperature is prescribed at \(T_d\) on \(\Gamma_D\)

Associated implementation:

`OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterThermal`

### Thermal analysis eigenvalues

The unknown temperature \(T\) of following thermal problem

can be computed using a spectral decomposition involving the following eigenvalue problem

where:

\(\mu_i\) is the i-th eigenvalue

\(T_i\) is the i-th eigenmode

\(T_d\) is the prescribed temperature

\(\rho\) is the material density

\(C_p\) is the thermal capacity

\(\lambda_c\) is the thermal conductivity

the temperature is prescribed at \(T_d\) on \(\partial \Omega\)

Associated implementation:

`OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterThermalEigenValues`

### Thermo linear elasticity analysis

We investigate a weakly coupled thermoelastic problem. As such, the first step is to perform a classic thermal analysis to find the temperature \(T\) (see thermal analysis for more details). The second step is, for a given temperature field, to find the displacement governed by the following problem

where:

\(\sigma^{elas}\) is the elastic contribution to the stress tensor, satisfying the Hooke’s law

\(e(u)\) is the strain tensor

\(\sigma^{ther}\) is the thermal contribution to the stress tensor, satisfying the following relation

\({\cal A}\) is the fourth order elasticity tensor given by

with \(\delta_{ij}\) the kronecker symbol.

\(\lambda\) and \(\mu\) are the Lamé coefficients

\(g\) is the density of surfaces forces

\(f\) is the density of volumic forces

the body is clamped on \(\Gamma_D\)

Associated implementation:

`OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterThermoElastic`

### Static buckling analysis

For this analysis, we refer to the following theoretical documentation 5.

### Modal analysis

We are looking for the eigenpulsations \(\omega_j > 0\) and its associated eigenmodes \(r_j\) such that

where:

\(\sigma\) is the stress tensor satisfying the Hooke’s law

\(\rho\) is material the density

\(j\) is the mode id

the body is clamped on \(\Gamma_D\)

Associated implementation:

`OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterModal`

### Harmonic analysis

We consider the harmonic response problem of a visco-elastic body depicted by the following problem

where:

\(\sigma\) is the stress tensor satisfying the Hooke’s law

\(\rho\) is the material density

\(c\) is the damping operator

\(g\) is the density of surfaces forces acting on \(\Gamma_N\)

the body is clamped on \(\Gamma_D\)

Associated implementation:

`OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterHarmonicForced`

### Physical analysis summary

Note that the following table only specify whether a specific analysis is available with OpenPisco when using a given physical solver. “False” does not mean it is not possible to actually run an analysis with a solver; at most it means the interface to do so is not available.

Analysis |
Code_Aster |
Freefem++ |
Zset |
Muscat solver |
---|---|---|---|---|

Linear elasticity analysis |
True |
True |
True |
True |

Thermo linear elasticity analysis |
True |
False |
False |
False |

Static buckling analysis |
True |
False |
False |
False |

Modal analysis |
True |
False |
False |
False |

Harmonic analysis |
True |
False |
False |
False |

Thermal analysis |
True |
False |
False |
False |

Thermal analysis eigenvalues |
True |
False |
False |
False |

## Fields available

### stress

We refer to the cauchy stress tensor, denoted by \(\boldsymbol{\sigma}\).

### strain

In the small deformation framework, the strain tensor is given by \(e(u)=(\nabla u + (\nabla u)^T)/2\) where \(u\) is a displacement.

### von_mises

where \(s_{ij}\) are the components of the Stress deviator tensor \(\boldsymbol{\sigma}^\text{dev}\) defined by

\(\boldsymbol{\sigma}^\text{dev} = \boldsymbol{\sigma} - \frac{\operatorname{tr}\left(\boldsymbol{\sigma}\right)}{3} \mathbf{I}\). where

### elastic_energy

The elastic energy is given by \(\frac{1}{2}\boldsymbol{\sigma}(u):e(u)\)

## Scalars available

### volume

The volume of the domain is given by \(\int_{\Omega}{ dx}\)

### mass

The volume of the domain is given by \(\int_{\Omega}{\rho dx}\) where \(\rho\) is the material density.

### int_elastic_energy

It is simply the elastic energy integrated on the whole domain, that is \(\int_{\Omega}{\frac{1}{2}\boldsymbol{\sigma}(u):e(u) dx}\)

### Fields/Scalars summary

Similarly to the physical analysis summary, an empty cell is to be understood as “the interface with this solver can not generate this field/scalar, it is not available as of now”.

We now denote by N,C,IP (ergo nodes, centroids, integration points) the field support type available for a field with respect to each solver. If any of this notation exist in a cell in the table below, it means the physical solver can generate this particular field on this specific support.

Analysis |
Code_Aster |
Freefem++ |
Zset |
Muscat |
---|---|---|---|---|

stress |
N,C,IP |
N |
||

strain |
N,C,IP |
|||

von_mises |
N,C,IP |
N |
||

elastic_energy |
N,C,IP |
N |
N |
N |

Next, we introduce the scalars in a similar way

Analysis |
Code_Aster |
Freefem++ |
Zset |
Muscat |
---|---|---|---|---|

volume |
False |
False |
False |
False |

mass |
False |
False |
False |
False |

int_elastic_energy |
True |
True |
True |
True |