Lagrange multipliers can be used to enforce constraints on the displacements of the system to machine precision. The introduction of these constraints change the structure and increase of the sparse matrix and restrict the type of solver to be used to only direct linear solvers.

It should be possible to implement this partitioned structure for displacement boundary conditions for non-linear problems and for using periodic boundary conditions for example. Tie constraints also naturally fall out of this.

A constraint matrix can be formulated for each of these describing the derivative of the gap function used for the particular constraint and provided under a different class of boundary conditions.

• Add interface class constraint
• Support constraint dofs in fem::mesh and fem::matrix
• Prohibit use of non-compatible linear solvers with zero-diagonal matrices

Implement handling of traction boundary conditions.

Plane strain is probably an easy example to get working and experiment with. A finite strain example will allow debugging of the three-dimensional example and set up internal variable store etc.

Develop the nonlinear geometry code

During long simulations it is useful to write out the state of the simulation such that it can be restarted from this checkpoint. The checkpoint should contain:

1. History variables (accumulated scalars / tensors)
2. Current displacements or nodal coordinates
3. Current time step and step size.

Required changes

The classes which hold state will require new constructors or logic to populate accumulation fields from restart data. Whether the restart should handle that different data as input is debatable. If there is a checkpoint at which point a simulation diverges due to bad material properties or time expiry on a shared compute device, it would be sensible to allow a change however this would break the continuity of the analysis and should perhaps be avoided entirely.

Each constitutive model should specify the internal variables which are listed as the history variables and provide this set through a constant method. The fem_submesh that is responsible for the internal variables is also responsible for serialisation of the scalar and tensor variables. Due to the potentially large size of the internal variable set, the results should be in binary format with the option of an ASCII output.

Classes with restart information could provide this in json format so this could be written / appended to a restart file <simulation_name>.restart with a mapping to a binary object containing the history variables. This is best handled in the fem*matrix classes as the time / load incrementation occurs here.

Serialisation libraries

The cereal library is header only and allows for a binary representation of the matrices. It seems somewhat easy to implement:

#include <Eigen/Dense>

#include <cereal/archives/binary.hpp>
#include <cereal/cereal.hpp>

#include <fstream>
#include <type_traits>

namespace cereal
{
template <class Archive, typename MatrixType>
inline std::enable_if_t<traits::is_output_serializable<BinaryData<typename MatrixType::Scalar>, Archive>::value> save(
    Archive& archive, MatrixType const& m)
{
    typename MatrixType::Index rows = m.rows();
    typename MatrixType::Index cols = m.cols();
    archive(rows);
    archive(cols);
    archive(binary_data(m.data(), rows * cols * sizeof(typename MatrixType::Scalar)));
}

template <class Archive, typename MatrixType>
inline std::enable_if_t<traits::is_input_serializable<BinaryData<typename MatrixType::Scalar>, Archive>::value> load(
    Archive& archive, MatrixType& m)
{
    typename MatrixType::Index rows, cols;
    archive(rows);
    archive(cols);

    m.resize(rows, cols);

    archive(binary_data(m.data(),
                        static_cast<std::size_t>(rows * cols * sizeof(typename MatrixType::Scalar))));
}
}

int main()
{
    Eigen::MatrixXd test1 = Eigen::MatrixXd::Random(10, 3);
    Eigen::MatrixXd test2 = Eigen::MatrixXd::Random(3, 5);

    {
        std::ofstream out("eigen.cereal", std::ios::binary);
        cereal::BinaryOutputArchive archive_o(out);
        archive_o(test1);
        archive_o(test2);
    }

    std::cout << "test1:" << std::endl << test1 << std::endl;
    std::cout << "test2:" << std::endl << test2 << std::endl;

    Eigen::MatrixXd test1_loaded, test2_loaded;
    {
        std::ifstream in("eigen.cereal", std::ios::binary);
        cereal::BinaryInputArchive archive_i(in);
        archive_i(test1_loaded, test2_loaded);
    }
    std::cout << "Norm of test matrix 1: " << (test1 - test1_loaded).norm() << "\n";
    std::cout << "Norm of test matrix 2: " << (test2 - test2_loaded).norm() << "\n";
}

and also supports containers in the standard library and json output (https://uscilab.github.io/cereal/index.html).

Implementation steps

• Have constitutive models populate a set of restart variables
• Print out current time information
• Create a checkpoint class to give to fem_submesh classes with file name information

Each fem_submesh owns and can print out their internal variables. This is because each submesh has a unique interpolation function and can printout based on this information alone.

Constitutive model

• Add a serialise() method using the code snippet.

This function would need to include its parent (information about the element or an index) and a method to construct its data with restart information, using pseudo c++:

class internal_variables
{
public:
    /// Serialise data to file
    void serialise(std::string const& file_tag,
                      std::set<variable::scalar> const& scalars
                      std::set<variable::second> const& tensors);

    /// Deserialise from file to memory
    void deserialise(std::string const& file_tag);
};

See 'A Framework for Finite Strain Elastoplasticity Based on Maximum Plastic Dissipation and the Multiplicative Decomposition, Part II'. This leads to a nice linear solve based on a diagonalised mass matrix rather than the average of the nodal values. Perhaps we can offer both methods?

If there is a numerical issue, the linear solver could compute a NaN and no error is currently thrown. To fix this wrap the linear solvers in floating point environment exceptions similar to the internal variable update.

I opened this issue is as a reminder to check.

There will be problems with integer overflow / undefined behaviour if there are more than roughly 2 billion elements or degrees of freedom. However the likelihood of such large meshes on a single node is low. We should check a large mesh and see if this will cause problems on a single node with high memory (~356 GB or there about) with iterative solvers.

An easy solution is to change everything 'scalable' to int64 and be done with it. This will require the linear solvers and associated ordering packages to be changed out to support 64 bit integers. The GPU solver will need changing because currently they are setup for 32 bit integers.

Implement the neo-Hookean weakly compressible form of the hyperelastic model.

Allow reaction forces at the nodes to be output. Aids in model sanity checking.

Tracker.

It should be possible to include the Gaussian chain in the microsphere model.

Develop code to perform continuous integration of the code.

This issue tracks the progress towards inclusion of MPI based parallelisation. Before this is included there should be a set of test codes to ensure that the currently used libraries can be enabled with MPI. Once this is completed, the linear solver backends need to be changed and the io implementation will be swapped out to include parallel file io.

For a worthwhile implementation it makes sense to restrict the choice of solvers to a 'fully' distributed such that no particular process has a complete view of the stiffness matrix. This method scales better because no one process has a world view on the problem.

1. Include external c++ wrappers for MPI (mpi-bind)
2. Distributed PaStiX
3. Distributed MUMPS (slight issues with this)
4. CPU Iterative solver suite
5. GPU Iterative solver suite
6. Distributed VTK IO

In order to avoid locking problems for incompressible models (such as the rubber models implemented), the weak form needs to be adjusted. This can be done through the Q1P0 family of elements, where the additional pressure variable is condensed out on the element level. However this is only valid for the discontinuous pressure field and higher order shape functions need to include the entire (u, p) function set and solve a partitioned system.

This is a tracking issue of the progress.

Implement a method to measure the code coverage. Aim for 50% code coverage for v0.1.

When the code is compiled with a debugging option, two tests fail GeometryTest and MatrixSystemTest.

Use boost::filesystem (or std::filesystem when stable) to output all the data into a visualisation directory. This avoids cluttering the mesh and simulation definition files

g++-7
libssl https
hwloc

Provide an option for fixed and adaptive time/load stepping for static problems. Fixed time stepping has only a queue and the fixed times. Non-convergence will result in a failed simulation.

Including multiphysics will be essential for interactions between mechanical and diffusion processes. This requires a robust way of transferring data between the modules to avoid excessive memory copying.

Requirements

This feature requires new infrastructure in the preprocessor to schedule the simulations. The requirements of this structure are:

• Propagation of results through multiple cases
• Apriori error checking. This will require construction of objects and checking inputs
• Arbitrary mixing of modules (mechanical followed by diffusion for example)
• Efficient access to internal variables and solution vectors to avoid excessive memory copies of large structures

Proposed implementation

An algorithm for scheduling the simulations is required, as a particular module may have multiple steps, where the steps reuse the previous simulation state. The scheduling model must be generic enough to handle a combination of procedures.

We need eigenvalue solvers in the mix. Add in the unsupported ARPACK wrapper. This library is used in Octave and other numerical packages so it should be robust.

Using structured bindings (decomposition bindings...?) and not using one of the variables results in lots of warnings. Using the [[maybe_unused]] specifier works in the latest clang and gcc. Once these new compilers hit the stable repos add the specifier to each call site where a variable is not used.

Implement test routines to test the functions and classes within the code.

Implement handling for displacement boundary conditions.

OpenMP looks like crap.

We should move towards c++11 standard threads or something a compatible with this library.

The main multithreaded portions of the program are inside the internal variables update and the matrix assembly. We can swap out those manual for loops for something like:

https://www.threadingbuildingblocks.org/tutorial-intel-tbb-generic-parallel-algorithms

Got some thoughts on this Shadi? We can still use an integer index which is handy when we have a lot of arrays to index into, and also means the code doesn't have to change too much.

There is a problem with the Tetrahedron4 element. The displacement field is correct but all the other fields like stress or damage are NAN. Tetrahedron10 works perfectly as you can see in the last picture.

Complete integration tests would be useful in ensuring the results are actually correct. This could be done using a python script that checks against analytical solutions.

Develop basic website for neon including documentation

Implement the small strain J2 plasticity.

The following library seems interesting and it can ease the effort put into JSON tests.
https://github.com/nlohmann/json

A trivial optimisation for the linear solver classes is to hold onto the state of the linear solver if the pattern hasn't changed. This was done for the LDLT PaStiX solver but it can be applied to all solvers.

This easy optimisation results in an impressive speed up of around 1/3 of the computation time on a unstructured grid with quadrature tetrahedrons, which is especially useful in non-linear solves where the sparse matrix pattern is unlikely to change.

Prism elements are the only missing common element from the library. We should add them in.

If we want to run fatigue computations, it will be easier to have a "cyclic" option that changes the interpolation of the boundary conditions to fit a sinusoidal wave. An input example may look like:

"Cyclic" : "True",
"Amplitude" : [0.12,0.10],
"Period" : [10,10],
"NumperOfCycles" : [2,1],
"Phase" : [0,0]

In this case "Times" and "Values" will be ignored and the respective values will be computed base on the given wave.

For surfaces of three-dimensional elements we need compute the mapping using the norm of the cross product.

Finite strain plasticity in three dimensions is kaputt. The reference by Neto provides an implementation for finite plasticity (J2) that we can compare against.

• Implement two dimensional theory and associated framework
• Test against the reference implementation
• Attempt the three dimensional theory again if it works

Fixed size numerical structures (vectors and matrices) are critical to extracting the maximum performance out of the code by providing as much information as possible to the optimiser. Depending on each simulation type, the matrices representing the material tangent, stress and deformation tensors are known. Some generic classes, for example the internal variables class, already take template parameters for the size of the rank two and rank four tensor. The dimension information (one, two or three) and the degrees of freedom can be encoded into these traits.

To deal with the quirks of each simulation type (axisymmetric, plate, shell, solid, beam, plane strain and plane stress) we should introduce a traits class that specifies this information for each simulation type. This generates a little additional work but should encapsulate the size information.

Firstly, we need to introduce a strongly-typed enum to provide a compile-time logic.

#include <type_traits>

enum class type : uint8_t {plane_stress, plane_strain, solid, plate, shell, axisymmetric};

template <type T, bool>
struct traits;

template <type::solid, bool is_symmetric = true>
struct traits
{
    static auto constexpr value_type = type::solid;
    static auto constexpr size = 3;

    using rank_two_tensor = matrix3;
    using rank_four_tensor = std::conditional<is_symmetric, matrix6, matrix9>::type;

    using internal_variables_type = internal_variables<rank_two_tensor, rank_four_tensor>;
    using constitutive_type = constitutive_model<type::solid, internal_variables_type>;
};

This should be able to handle the unsymmetrical 9x9 constitutive models for three-dimensional theory if required through the template parameter. In addition, if analysis specific post-processing steps are required, further template specialisations can be used for the parts which are not common.

A downside of this approach is the additional boiler plate and the fact that most of the implemented code will still be dependent on the type with what code can be reused, will be.

Probably should get this spun up.

Pressure loads follow the deformation of the body. This introduces another part of the stiffness matrix called the 'load stiffness' which causes the system to lose its symmetry. Issue #11 allows the linear solver to select an appropriate solver if the system is unsymmetric.

• Add follower_pressure class with external_stiffness and external_force interface

This is to track the implementation of the crack closure effect in the Chaboche damage model

For some exotic material models (i.e. this rubber ageing) there are some problems with defining a sensible tolerance. We should add the option in for using an absolute or relative tolerance from the input file.

The routines where the mesh is gathered and written to VTK format have multiple distinct tasks. It should be possible to spawn asynchronous tasks (with std::async) to fill the mesh data structure, extrapolate nodal variables and add these into the VTK container with the appropriate locks. This would ensure that file io is loosely parallelised and should net a speedup for larger models.

This asynchronous approach could also be used in other parts of the code (preprocessing)

Add citations to the doxygen so model origins can be checked / traced.

Develop the nonlinear material code

Displacement controlled loadings as they currently stand cause spurious plastic deformations in the elements connected to displacement nodes. This is because the load is applied such that high localised strains are computed at the quadrature points.

It should be that a global or local linear solution is computed first as an approximation to the deformation state.

The update of the internal variables with a non-zero time step is done before invoking perform_equilibrium_iterations and it is done again inside perform_equilibrium_iterations with a zero time step. However, for visco-plastic materials, the update with a non-zero time step should be done in the first N-R iteration only.

The affine microsphere model now includes ageing effects but for purely mechanical studies this is overkill. Split these classes and include ageing in a derived class

Currently this is not enforced and needs to occur at time of construction of the fem*Matrix. The mesh objects should be able to report if they will produce an unsymmetric matrix based on their constitutive model / boundary conditions.

Currently the input file specifies the choice of linear solver as "ConjugateGradient" or "SparseLU", but this choice can conflict with the choice of constitutive model, constraint enforcement and line search methods. A better approach would be to specify the class of solvers and have the code automatically deduce the best approach based on:

• Constitutive model
• Nonlinear solver options
• Constraint equations
• Boundary conditions (follower loads)

The suggested supported solvers options could be

• "PaStiX" (LU, LLT) direct solver
• "MUMPS" (LU, LLT) direct solver
• "Direct" (Eigen inbuilt LU, LLT)
• "Iterative" (Eigen CG and BiCGStab (or GMRES))
• "IterativeGPU" (CG and BiCGStab on the GPU)
• "DirectGPU" (QR factorisation from CUDA library)

Changes required to be implemented are:

• Mesh classes query their objects to determine symmetry status
• femMatrix classes check their nonlinear options and meshes and deduce correct solver

The class structure of linear solvers should follow the corresponding hierarchy with two classes Iterative and Direct serving as base classes with the derived classes containing specialisations (iterative have tolerances and direct solvers with saving of the factorised matrices). For iterative solvers, there are options to specify the precondition to improve convergence properties. These could also be specified when possible in a solver agnostic way, with "Preconditioner" : "Diagonal" or "Preconditioner" : "IncompleteFactorisation", or with a similar method.

An implementation of these changes should increase the linear solver efficiency and avoid incorrect results.

WIth c++17 a variant type (a type safe union) has been introduced into the standard library. Instead of using inheritance and polymorphic pointers, the boundary conditions can be modelled as a union of types with an .external_force() parameter and accessed through a visitor and a lambda.

For this, we need to split the boundary conditions into two types; those who contribute to the right hand side and those who contribute terms in the stiffness matrix and terms to the right hand side.

The variant type can be defined as

using boundary_type = std::variant<Traction, Pressure, BodyForce>;

A factory method can allocate the correct type to the variant

// Construct a list of boundary conditions
std::vector<boundary_type> boundary_conditions;
boundary_conditions.emplace_back(Pressure{});
boundary_conditions.emplace_back(Traction{});
boundary_conditions.emplace_back(BodyForce{});

Within the system matrix assembler, the original code requires a small modification

for (auto const& boundary_condition : boundary_conditions)
{
    std::visit([&](auto const& boundary) {
        // Assemble boundary condition contribution for the right hand side
        // Loop over elements
        for (auto element = 0; element < elements; ++element)
        {
        }
    }, boundary_condition);
}

From the input file, the traction boundary conditions are usually given as time history for each degree of freedom:

"x" : [1.0, 2.0, 3.0],
"y" : [0.0, 1.0, 2.0],
"z" : [0.0, 0.0, 0.0]

This format needs to be stored internally for operation. Currently, each traction boundary is held in a std::vector inside a container class. An advantage is that each degree of freedom for the traction is handled separately and does not contain a notion of what degree of freedom it represents. This abstracts all code and allows greater flexibility while simplifying the base class.

A downside is that this approach causes dynamic memory allocation. It would be easier to have a std::array<> and set a boolean flag for active or inactive DoFs.

For handling the boundary conditions the following framework is proposed. We are primarily concerned with a few common types for non-follower loads to begin with

• Traction
• Body force (gravity)
• Pressure

Traction boundary condition are composed of three DoFs in the most general case. Between time/load steps, the traction value needs to be interpolated. This can be generalised to a Vector3 type.

Each boundary condition needs to have an encoding between the given DoF, "x" for example, and the corresponding numerical number. This data is specific to the physics being simulated. A solid mechanics simulation may have a different ordering to an axisymmetric model and so forth.