dbeurle / neon Goto Github PK
View Code? Open in Web Editor NEWA finite element code
License: Other
A finite element code
License: Other
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.
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.
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
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?
Implement handling for displacement boundary conditions.
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.
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.
constraint
fem::mesh
and fem::matrix
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.
The following library seems interesting and it can ease the effort put into JSON tests.
https://github.com/nlohmann/json
Allow reaction forces at the nodes to be output. Aids in model sanity checking.
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:
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.
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).
fem_submesh
classes with file name informationEach 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.
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);
};
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.
Implement a method to measure the code coverage. Aim for 50% code coverage for v0.1.
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.
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.
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.
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 neo-Hookean weakly compressible form of the hyperelastic model.
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.
Develop the nonlinear geometry code
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.
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.
Implement test routines to test the functions and classes within the code.
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.
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.
For surfaces of three-dimensional elements we need compute the mapping using the norm of the cross product.
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.
follower_pressure
class with external_stiffness
and external_force
interfaceThis is to track the implementation of the crack closure effect in the Chaboche damage model
Prism elements are the only missing common element from the library. We should add them in.
Probably should get this spun up.
Implement the small strain J2 plasticity.
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)
Finite strain plasticity in three dimensions is kaputt. The reference by Neto provides an implementation for finite plasticity (J2) that we can compare against.
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.
Implement handling of traction boundary conditions.
Develop code to perform continuous integration of the code.
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 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.
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.
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
Add citations to the doxygen so model origins can be checked / traced.
When the code is compiled with a debugging option, two tests fail GeometryTest
and MatrixSystemTest
.
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:
The suggested supported solvers options could be
Changes required to be implemented are:
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.
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.
g++-7
libssl https
hwloc
Tracker.
It should be possible to include the Gaussian chain in the microsphere model.
Develop the nonlinear material code
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.
This feature requires new infrastructure in the preprocessor to schedule the simulations. The requirements of this structure are:
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.