SimVascular / svFSIplus

svFSIplus is an open-source, parallel, finite element multi-physics solver.
https://simvascular.github.io/documentation/svfsi.html
Other
10 stars 26 forks source link

Object-oriented solid material models #178

Open mrp089 opened 9 months ago

mrp089 commented 9 months ago

Use Case

This affects all solid material models that are currently implemented and future ones.

Problem

We currently have a massive switch/case statement for our solid materials, inherited from svFSI. Many of these material models overlap heavily. Furthermore, we cannot combine material models to form a new one, making it hard to implement new (more complex) material models (like G&R).

Solution

Create a class structure for all existing material models.

Alternatives considered

FEBio has a comprehensive collection of material models. This can help us figure out our material models and some of their inheritance structures.

Abaqus has interfaces for user material subroutines. There is one for general solid materials and an invariant-based one for hyperelastic materials.. This can help to define interfaces between abstract (e.g., hyperelastic) and derived classes (e.g., Neo-Hooke, Mooney-Rivlin).

Additional context

This is a good opportunity to take "inventory" of what materials are currently implemented. We can use our (currently under-utilized) Doxygen documentation when implementing new classes to document the governing equations. This also helps visualize the inheritance structure (see here a 0D example).

Code of Conduct

lpapamanolis commented 9 months ago

I have started looking into this: I think the hierarchical structure implemented in Abaqus is very close to how we would like to restructure our code. The rough class hierarchy would look something like this:

  1. Material: user-defined label, user should be able to combine different pre-defined material properties/models from the levels below
  2. (Material) Properties: Material properties grouped in broad categories (categories in abaqus: elastic mechanical properties, inelastic mechanical properties, thermal properties, acoustic properties, hydrostatic fluid properties, equations of state, mass diffusion properties, electrical properties). Here we could add "growth and remodelling properties".
  3. (Material) Model Class: Classes of models that inherit properties from level 2. For example, within the elastic mechanical properties category we can have linear elastic, viscoelastic, poroelastic, hyperelastic models etc.
  4. (Material) Model: Specific model used from a given property/class from levels 2/3. For example, hyperelastic models can include St.Venant-Kirchhoff, modified St.Venant-Kirchhoff, Neo-Hookean, Mooney-Rivlin, HGO (Holzapfel-Gasser-Ogden), Guccione, and Holzapfel-Ogden models.
  5. (Model) Parameters: Specific numerical values for the parameters of the chosen models, or other model-specific variables.

For example, suppose we want to model an arterial wall (domain 1) externally supported by a sheath (domain 2). Then we can define material_wall=elastic.hyperelastic.neo_hookean, and material_sheath=elastic.linear.standard for the sheath. We can then define parameters material_wall.shear_modulus and material_wall.bulk_modulus for the Neo_hookean model and so forth for the linear elastic model.

If we want to add growth to the arterial wall, then we would just combine it in the same variable: material_wall=gnr.costrained_mixture.equilibrated and add parameters material_wall.add_constituent(constituent_parameters).

Of course not everything can be combined with everything else, but Abaqus already has documentation on what can be combined with what, under which conditions, which building blocks are incompatible etc. For growth and remodeling, we also need to think of an appropriate hierarchy based on the above template, and if/how it should be combined with other models.

I am currently looking into and documenting what models are already implemented in svFSI+ and will follow-up on that. Let me know if you have any comments on the above in the meantime.

ktbolt commented 9 months ago

@lpapamanolis The hierarchical structure you propose looks like a nice way to identify material model components that can be mapped to classes.

An important implementation goal to keep in mind is enabling code reuse. This can be implemented using inheritance or composition. composition is typically preferred because it leads to less coupled classes and is easier to maintain.

lpapamanolis commented 9 months ago

I'm leaving here a rough schematic of the structure I proposed above for future reference. We can iterate and flesh it out along the way.

image
yuecheng-yu commented 7 months ago

Thanks for discussion with @mrp089 and @MatteoSalvador yesterday. We observed that to implement Object-Oriented Solid Material Model and further to implement G&R Constraint Mixture Constituent Model, we may need material point to store the material information at every gauss point. The purpose is similar to what currently class stModelType does in get_pk2cc(): load the material type and properties. However, class stModelType is inherited from class dmnType: com_mod.eq[cEq].dm[cDmn].stM.

As a start, I will write a class MaterialPoint and create "material points" at each gauss point. I think it should be inherited from class mshType: com_mod.msh[cMsh]. During initialization, "material points" will copy the material related information from the domain-based class stModelType. At each Newton Iteration, "material points" will manage all solid material behaviors, such as store/retrieve/update material properties, compute stress/tangent at each gauss point.

One concern about the above implementation I can imagine is the efficiency and MPI compatibility. Because the class MaterialPoint requires to create an instance at every gauss point using heap memory, I am not sure if it will slow the program or have any issue with MPI. It should be firstly tested after implementation.

@ktbolt @lpapamanolis Any comments or thoughts are very welcomed!

ktbolt commented 7 months ago

@yuecheng-yu A good start! Your material model work should probably be in a separate Issue.

A few comments.

Saying that stModelType class is inherited from the dmnType class means

class stModelType : public dmnType {}

The dmnType is a composition of several classes of which stModelType is one.

Before you write a single line of code you should create a design document that defines what a MaterialPoint object is (what it is used for, what its data looks like, etc), what it does and how it interacts with other parts of the code (API). This document will be part of the code design documents available to developers. Have a look at the svZeroDSolver documentation and see how the code architecture and implementation is almost totally opaque to a new developer because of no design documents.

You can then create an implementation design that sketches out how the MaterialPoint object will be implemented. You can do this using pseudo code. This is a good place to understand data layout and how memory will be used, will address any issue with memory usage (e.g. do you need to write data to disk).

You can then prototype your implementation in a branch that the team can use for code reviews.

The MaterialPoint data will need to be partitioned across processors like other data (e.g. mesh).