new loads

️ More versatile fixation system

It is now possible to define more complex sets of fixations.

Previously, 3 sets of nodes were hard-coded in the input file as "contact_pts", "axis_pt1", "axis_pt2".

Example:

    p['contact_pts'] = [[-10.20362, -17.46838, -229.9061],
                        [-11.92466, 26.3042, -229.5354]]
    p['axis_pt1'] = f'{path}/LTMJ.off'
    p['axis_pt2'] = [-8.716309, 79.13171, 233.8385]
    p['fixations'] = {
        'contact_pts': ['x'],
        'axis_pt1': ['x', 'y', 'z'],
        'axis_pt2': ['x', 'z']
    }

In this new version, the user can define any number of fixation sets and name them as he/she wants. It is thus possible to define much more complex senarios than before.

The previous commands should be gathered in a single list of fixations:

    p['fixations'] = [
        {
            'name': 'contact_pts',
            'nodes': [[-10.20362, -17.46838, -229.9061],
                        [-11.92466, 26.3042, -229.5354]],
            'direction': 'x'
        },
        {
            'name': 'axis_pt1',
            'nodes': f'{path}/LTMJ.off',
            'direction': ['x', 'y', 'z']
        },
        {
            'name': 'axis_pt2',
            'nodes': [-8.716309, 79.13171, 233.8385],
            'direction': ['x', 'z']
        }
    ]

This list of fixations can be of any size.

️ In consequence, the old syntax of input files is not compatible with this new version of fossils! (but translation is rather straightforward).

️ Additional loads

Additional loads can be added with a syntax similar to the new syntax for fixations:

    p['loads'] = [
        {
            'name': 'force1',
            'nodes': [-10.20362, -17.46838, -229.9061],
            'values': [0., 0., -100.]
        },
        {
            'name': 'force2',
            'nodes': [-11.92466, 26.3042, -229.5354],
            'values': [0., 0., 100.]
        }      
    ]

with:

• name: a user-defined string, mainly used to name the group of nodes,
• nodes: as for fixations, it can be a point defined by its 3 coordinates between square brackets ([x, y, z]), or a list of points defined in the same manner (e.g. [[x1, y1, z1], [x2, y2, z2], ...]) or a mesh file containing one or more vertices.
• values: the force vector applied to each node.

A new model named models/dolicorhynchops/dolicorhynchops_10k_extloads.py has been created to test this new feature.

️ New load distribution models on surface meshes

In addition to the "Ad Hoc Uniform Traction Model" ('U'), the "Tangential-Traction Model" (T) and the "Tangential-Plus-Normal-Traction Model" (T+N), the user can now use 2 new load distributions:

Directional-Traction Model (D)

The forces are distributed on the nodes of a given surface (file key). They are all pointing in the focal_pt direction seen as a vector. They are scaled so that the total resulting force has the magnitude prescribed by the user (with the key force)

Example:

    p['muscles'] = [
        {
            'method': 'D',
            'file': f'{path}/Lmuscle.stl',  # surface mesh
            'force': 100.,                  # total force
            'focalpt': [-1, 1, 1],          # <= used as a direction
        }
    ]

Note: the keyword focal_pt can be misleading here. The vector given as argument is a direction. However we have kept the string focal_pt as a keyword to have the same keywords for all the methods.

:dart:Normal-Traction Model (N)

In this new method, the nodal forces are applied in the perpendicular direction to each facet of the surface mesh. The focalpt keyword is ignored. Forces are scaled to have the prescribed resulting magnitude (the sum of all the forces is force).

Example:

    p['muscles'] = [
        {
            'method': 'N',
            'file': f'{path}/Lmuscle.stl',  # surface mesh
            'force': 100.,                  # total force
            'focalpt': [0, 0, 0],           # <= ignored
        }
    ]

cxxfem/src/femSolver.cpp

@@ -67,8 +67,8 @@ void Solver::solve()
 
     create_partition();
     fill_partition();
-    fill_rhs();
     fill_lhs();
+    fill_rhs();
     apply_bc();
     solve_static();
 
@@ -206,11 +206,11 @@ void Solver::fill_partition()
     std::cout << "\telasped = " << timers["fill_partition"].read() << "s\n";
 }
 
-/// fill right-hand side of the system of equations
+/// fill left-hand side of the system of equations
 
-void Solver::fill_rhs()
+void Solver::fill_lhs()
 {
-    timers["fill_rhs"].start();
+    timers["fill_lhs"].start();
 
     std::cout << "assembly of LHS...\n";
 
@@ -527,15 +527,15 @@ void Solver::fill_rhs()
         for (size_t j = 0; j < Kijv[i].size(); ++j)
             delete Kijv[i][j];
 
-    timers["fill_rhs"].stop();
-    std::cout << "\telasped = " << timers["fill_rhs"].read() << "s\n";
+    timers["fill_lhs"].stop();
+    std::cout << "\telasped = " << timers["fill_lhs"].read() << "s\n";
 }
 
 /// compute right-hand side 'f' (sources)
 
-void Solver::fill_lhs()
+void Solver::fill_rhs()
 {
-    timers["fill_lhs"].start();
+    timers["fill_rhs"].start();
 
     std::cout << "assembly of RHS...\n";
 
@@ -627,8 +627,8 @@ void Solver::fill_lhs()
         }
     }
 
-    timers["fill_lhs"].stop();
-    std::cout << "\telasped = " << timers["fill_lhs"].read() << "s\n";
+    timers["fill_rhs"].stop();
+    std::cout << "\telasped = " << timers["fill_rhs"].read() << "s\n";
 
 }