Up elastic collision algo

Implement contact model for elastic bodies #2203

feelpp/quickstart/CMakeLists.txt

@@ -78,6 +78,7 @@ feelpp_add_application( elasticity_2d SRCS qs_elasticity.cpp INCLUDE_IN_ALL DEFS
 feelpp_add_application( elasticity_3d SRCS qs_elasticity.cpp INCLUDE_IN_ALL DEFS FEELPP_DIM=3 EXEC QS_ELASTICITY_3D INSTALL )
 feelpp_add_application( elasticity_pure_traction_2d SRCS qs_elasticity_pure_traction.cpp INCLUDE_IN_ALL DEFS FEELPP_DIM=2 EXEC QS_ELASTICITY_PURE_TRACTION_2D INSTALL )
 feelpp_add_application( elasticity_pure_traction_3d SRCS qs_elasticity_pure_traction.cpp INCLUDE_IN_ALL DEFS FEELPP_DIM=3 EXEC QS_ELASTICITY_PURE_TRACTION_3D INSTALL  )
+feelpp_add_application( elasticity_contact SRCS qs_elasticity_contact.cpp INCLUDE_IN_ALL DEFS EXEC elasticity_contact INSTALL )
 
 feelpp_add_testcase( elasticity PREFIX qs CATEGORY quickstart )
 add_dependencies(quickstart ${QS_ELASTICITY_2D} ${QS_ELASTICITY_3D})

feelpp/quickstart/elasticity/contact/ball.cfg

@@ -0,0 +1,11 @@
+directory=qs_elasticity_contact/ball
+specs=$cfgdir/ball.json
+
+pc-type = lu
+
+[gmsh]
+geo-variables-list=dim=2
+
+[exporter]
+static=true
+geometry=static

feelpp/quickstart/elasticity/contact/ball.geo

@@ -0,0 +1,13 @@
+h = 0.8;
+Point(1) = {0, 0, 0, h};
+Point(2) = {20, 0, 0, h};
+Point(3) = {-20, 0, 0, h};
+
+Circle(1) = {3, 1, 2};
+Circle(2) = {2, 1, 3};
+Curve Loop(1) = {1, 2};
+Plane Surface(1) = {1};
+
+
+Physical Curve("Wall", 4) = {1,2};
+Physical Surface("Caoutchouc", 5) = {1}; 

feelpp/quickstart/elasticity/contact/ball.json

@@ -0,0 +1,111 @@
+{
+    "Name": "2D falling ball",
+    "ShortName": "FallingBall2D",
+    "Models": {
+        "LinearElasticity": {
+            "Name": "FallingBall2D",
+            "dimension": 2,
+            "order": 2,
+            "Materials":[
+                "Caoutchouc"
+            ],
+            "loading": {
+                "F1": {
+                    "type" : "Gravity",
+                    "parameters" : {
+                        "expr": "{0,-0.1}",
+                        "unit": "N"
+                    }
+                }
+            }
+        }
+    },
+    "Materials": {
+        "Caoutchouc": {
+            "type": "LinearElastic",
+            "parameters": {
+                "E": {
+                    "value":"75.0",
+                    "unit":"Pa",
+                    "description":"Young's modulus of caoutchouc",
+                    "min": "75.0",
+                    "max": "75.0"
+                },
+                "nu": {
+                    "value":"0.25",
+                    "unit":"",
+                    "description":"Poisson ratio of caoutchouc",
+                    "min": "0.25",
+                    "max": "0.25"
+                },
+                "rho": {
+                    "value":"1",
+                    "unit":"kg/m^3",
+                    "description":"Density of caoutchouc",
+                    "min": "1",
+                    "max": "1"
+                }
+            }
+        }
+    },
+    "Meshes": {
+        "LinearElasticity": {
+            "Import": {
+                "filename": "$cfgdir/ball.geo",
+                "partition": 0,
+                "h": 0.4
+            }
+        }
+    },
+    "TimeStepping":
+    {
+        "LinearElasticity" :{
+            "steady": false,
+            "order" : 2,
+            "start": 0.0,
+            "end": 15.0,
+            "step": "1/10"
+        }
+    },
+    "InitialConditions": {
+        "LinearElasticity": {
+            "displacement": {
+                "expr": "{0.0,0.0}"
+            }
+        }
+    },
+    "BoundaryConditions": {
+        "LinearElasticity": {
+
+        }
+    },
+    "Collision":{
+        "LinearElasticity":{
+            "distance":-24,
+            "direction":"{0,-1}",
+            "theta":1,
+            "gamma0":200,
+            "fixedPoint":1,
+            "tolerance":0.001
+        }
+    },
+    "Parameters": {
+        "LinearElasticity": {
+            "E": {
+                "expr": 75
+            },
+            "nu": {
+                "expr": 0.25
+            },
+            "rho": {
+                "expr": 1.0
+            },
+            "beta": {
+                "expr": 0.25
+            },
+            "gamma": {
+                "expr": 0.5
+            }
+        }
+    }
+}

feelpp/quickstart/qs_elasticity_contact.cpp

@@ -0,0 +1,54 @@
+#include "qs_elasticity_contact.hpp"
+
+int main(int argc, char** argv)
+{
+    using namespace Feel;
+    int status;
+
+    try
+    {
+        Environment env(_argc = argc, _argv = argv,
+                        _desc = makeOptions(),
+                        _about = about(_name = fmt::format("elasticity_contact"),
+                                       _author = "Feel++ Consortium",
+                                       _email = "feelpp@cemosis.fr"));
+
+        auto jsonfile = removeComments(readFromFile(Environment::expand(soption("specs"))));
+        std::istringstream istr(jsonfile);
+        json specs = json::parse(istr);
+
+        if ( specs["/Models/LinearElasticity"_json_pointer].empty() )
+            throw std::runtime_error("No LinearElasticity model specified in the input file");
+
+        int dimension = specs["/Models/LinearElasticity/dimension"_json_pointer];
+        int order = specs["/Models/LinearElasticity/order"_json_pointer];
+
+        if ( dimension == 2 && order == 1)
+        {
+            ElasticContact<2, 1> ElasticContact(specs);
+            ElasticContact.run();
+        }
+        else if ( dimension == 2 && order == 2)
+        {
+            ElasticContact<2, 2> ElasticContact(specs);
+            ElasticContact.run();
+        }
+        else if ( dimension == 3 && order == 1)
+        {
+            ElasticContact<3, 1> ElasticContact(specs);
+            ElasticContact.run();
+        }
+        else if ( dimension == 3 && order == 2)
+        {
+            ElasticContact<3, 2> ElasticContact(specs);
+            ElasticContact.run();
+        }
+        else
+            throw std::runtime_error(fmt::format("Invalid dimension {} specified in the input file", dimension));
+    }
+    catch (...)
+    {
+        handleExceptions();
+    }
+    return 0;
+}