1 /// This is a simple numerical example of modeling a muscle fiber using the nonlinear elasticity solver
2 /// in a mixed displacement-pressure formulation. It roughly corresponds to Example 5.1 from the following paper:
3 /// M.H.Gfrerer and B.Simeon "Fiber-based modeling and simulation of skeletal muscles" 2020
4 ///
5 /// Author: A.Shamanskiy (2016 - ...., TU Kaiserslautern)
6 #include <gismo.h>
7 #include <gsElasticity/gsElasticityAssembler.h>
8 #include <gsElasticity/gsIterative.h>
9 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
10
11 using namespace gismo;
12
main(int argc,char * argv[])13 int main(int argc, char* argv[]){
14
15 gsInfo << "This is a muscle fiber benchmark with a mixed nonlinear elasticity solver.\n";
16
17 //=====================================//
18 // Input //
19 //=====================================//
20
21 std::string filename = ELAST_DATA_DIR"/muscleBeamMP.xml";
22 real_t youngsModulus = 3.0e5; // shear modulus 1e5;
23 real_t poissonsRatio = 0.5;
24 real_t density = 9e2;
25 real_t gravityAcc = -9.8;
26 // spatial discretization
27 index_t numUniRefDirX = 4;
28 index_t numUniRef = 0;
29 index_t numDegElev = 0;
30 bool subgridOrTaylorHood = false;
31 // output
32 index_t numPlotPoints = 1000;
33
34 // minimalistic user interface for terminal
35 gsCmdLine cmd("This is a muscle fiber benchmark with mixed nonlinear elasticity solver.");
36 cmd.addInt("x","xrefine","Number of uniform refinement along the beam axis",numUniRefDirX);
37 cmd.addInt("r","refine","Number of uniform refinement applications",numUniRef);
38 cmd.addInt("d","degelev","Number of degree elevation applications",numDegElev);
39 cmd.addSwitch("e","element","True - subgrid, false - TH",subgridOrTaylorHood);
40 cmd.addInt("s","points","Number of points to plot to Paraview",numPlotPoints);
41 try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
42 gsInfo << "Using " << (subgridOrTaylorHood ? "Taylor-Hood " : "subgrid ") << "mixed elements.\n";
43
44 //=============================================//
45 // Scanning geometry and creating bases //
46 //=============================================//
47
48 // scanning geometry
49 gsMultiPatch<> geometry;
50 gsReadFile<>(filename, geometry);
51
52 // creating bases
53 gsMultiBasis<> basisDisplacement(geometry);
54 gsMultiBasis<> basisPressure(geometry);
55 for (index_t i = 0; i < numDegElev; ++i)
56 {
57 basisDisplacement.degreeElevate();
58 basisPressure.degreeElevate();
59 }
60 for (index_t i = 0; i < numUniRef; ++i)
61 {
62 basisDisplacement.uniformRefine();
63 basisPressure.uniformRefine();
64 }
65 for (size_t p = 0; p < geometry.nPatches(); ++p)
66 for (index_t i = 0; i < numUniRefDirX; ++i)
67 {
68 static_cast<gsTensorNurbsBasis<3,real_t> &>(basisDisplacement.basis(p)).knots(0).uniformRefine();
69 static_cast<gsTensorNurbsBasis<3,real_t> &>(basisPressure.basis(p )).knots(0).uniformRefine();
70 }
71 // additional displacement refinement for stable mixed FEM
72 if (!subgridOrTaylorHood) // subgrid
73 basisDisplacement.uniformRefine();
74 else // Taylor-Hood
75 basisDisplacement.degreeElevate();
76
77 //=============================================//
78 // Setting loads and boundary conditions //
79 //=============================================//
80
81 // boundary conditions
82 gsBoundaryConditions<> bcInfo;
83 for (size_t p = 0; p < geometry.nPatches(); ++p)
84 for (index_t d = 0; d < 3; ++d)
85 {
86 bcInfo.addCondition(p,boundary::west,condition_type::dirichlet,nullptr,d);
87 bcInfo.addCondition(p,boundary::east,condition_type::dirichlet,nullptr,d);
88 }
89 // source function, rhs
90 gsConstantFunction<> gravity(0.,0.,gravityAcc*density,3);
91
92 //=============================================//
93 // Solving //
94 //=============================================//
95
96 // creating assembler for the displacement-pressure formulation
97 gsElasticityAssembler<real_t> assembler(geometry,basisDisplacement,basisPressure,bcInfo,gravity);
98 assembler.options().setReal("YoungsModulus",youngsModulus);
99 assembler.options().setReal("PoissonsRatio",poissonsRatio);
100 gsInfo << "Initialized system with " << assembler.numDofs() << " dofs.\n";
101
102 // setting Newton's method
103 gsIterative<real_t> solver(assembler);
104 solver.options().setInt("Verbosity",solver_verbosity::all);
105 solver.options().setInt("Solver",linear_solver::LDLT);
106
107 gsInfo << "Solving...\n";
108 gsStopwatch clock;
109 clock.restart();
110 solver.solve();
111 gsInfo << "Solved the system in " << clock.stop() <<"s.\n";
112
113 //=============================================//
114 // Output //
115 //=============================================//
116
117 // displacement and pressure as isogeometric fields
118 gsMultiPatch<> displacement,pressure;
119 assembler.constructSolution(solver.solution(),solver.allFixedDofs(),displacement,pressure);
120 // construct stress field
121 gsPiecewiseFunction<> stresses;
122 assembler.constructCauchyStresses(displacement,pressure,stresses,stress_components::von_mises);
123
124 if (numPlotPoints > 0) // visualization
125 {
126 // constructing an IGA field (geometry + solution)
127 gsField<> displacementField(geometry,displacement);
128 gsField<> pressureField(geometry,pressure);
129 gsField<> stressField(assembler.patches(),stresses,true);
130 // creating a container to plot all fields to one Paraview file
131 std::map<std::string,const gsField<> *> fields;
132 fields["Displacement"] = &displacementField;
133 fields["Pressure"] = &pressureField;
134 fields["von Mises"] = &stressField;
135 gsWriteParaviewMultiPhysics(fields,"muscleBeam",numPlotPoints);
136 gsInfo << "Open \"muscleBeam.pvd\" in Paraview for visualization.\n";
137 }
138
139 return 0;
140 }
141