1 /// This is the structural solver benchmark CSM3 from this paper:
2 /// "Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow"
3 /// Stefan Turek and Jaroslav Hron, <Fluid-Structure Interaction>, 2006
4 ///
5 /// Author: A.Shamanskiy (2016 - ...., TU Kaiserslautern)
6 #include <gismo.h>
7 #include <gsElasticity/gsElasticityAssembler.h>
8 #include <gsElasticity/gsMassAssembler.h>
9 #include <gsElasticity/gsElTimeIntegrator.h>
10 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
11
12 using namespace gismo;
13
writeLog(std::ofstream & file,const gsMultiPatch<> & displacement,real_t simTime,real_t compTime,real_t numIters)14 void writeLog(std::ofstream & file, const gsMultiPatch<> & displacement,
15 real_t simTime, real_t compTime, real_t numIters)
16 {
17 // evaluating displacement at the point A
18 gsMatrix<> point(2,1);
19 point << 1., 0.5;
20 gsMatrix<> dispA = displacement.patch(0).eval(point);
21
22 // print: simTime dispAx dispAy compTime numIters
23 file << simTime << " " << dispA.at(0) << " " << dispA.at(1) << " "
24 << compTime << " " << numIters << std::endl;
25 }
26
main(int argc,char * argv[])27 int main(int argc, char* argv[]){
28 gsInfo << "Benchmark CSM3: dynamic deflection of an elastic beam.\n";
29
30 //=====================================//
31 // Input //
32 //=====================================//
33
34 std::string filename = ELAST_DATA_DIR"/flappingBeam_beam.xml";
35 real_t poissonsRatio = 0.4;
36 real_t youngsModulus = 1.4e6;
37 real_t density = 1.0e3;
38 real_t loading = 2.;
39 // space discretization
40 index_t numUniRef = 3;
41 index_t numDegElev = 0;
42 // time integration
43 real_t timeSpan = 2;
44 real_t timeStep = 0.01;
45 // output
46 index_t numPlotPoints = 1000;
47
48 // minimalistic user interface for terminal
49 gsCmdLine cmd("Benchmark CSM3: dynamic deflection of an elastic beam.");
50 cmd.addReal("l","load","Gravitational loading acting on the beam",loading);
51 cmd.addInt("r","refine","Number of uniform refinement application",numUniRef);
52 cmd.addInt("d","degelev","Number of degree elevation application",numDegElev);
53 cmd.addReal("t","time","Time span, sec",timeSpan);
54 cmd.addReal("s","step","Time step, sec",timeStep);
55 cmd.addInt("p","points","Number of sampling points to plot to Paraview",numPlotPoints);
56 try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
57
58 //=============================================//
59 // Scanning geometry and creating bases //
60 //=============================================//
61
62 // scanning geometry
63 gsMultiPatch<> geometry;
64 gsReadFile<>(filename, geometry);
65
66 // creating bases
67 gsMultiBasis<> basisDisplacement(geometry);
68 for (index_t i = 0; i < numDegElev; ++i)
69 basisDisplacement.degreeElevate();
70 for (index_t i = 0; i < numUniRef; ++i)
71 basisDisplacement.uniformRefine();
72
73 //=============================================//
74 // Setting loads and boundary conditions //
75 //=============================================//
76
77 // boundary conditions
78 gsBoundaryConditions<> bcInfo; // numbers are: patch, function pointer (nullptr) for displacement, displacement component
79 bcInfo.addCondition(0,boundary::west,condition_type::dirichlet,0,0);
80 bcInfo.addCondition(0,boundary::west,condition_type::dirichlet,0,1);
81
82 // gravity, rhs
83 gsConstantFunction<> gravity(0.,loading*density,2);
84
85 //=============================================//
86 // Setting assemblers and solvers //
87 //=============================================//
88
89 // creating stiffness assembler
90 gsElasticityAssembler<real_t> assembler(geometry,basisDisplacement,bcInfo,gravity);
91 assembler.options().setReal("YoungsModulus",youngsModulus);
92 assembler.options().setReal("PoissonsRatio",poissonsRatio);
93 assembler.options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff);
94 gsInfo << "Initialized system with " << assembler.numDofs() << " dofs.\n";
95
96 // creating mass assembler
97 gsMassAssembler<real_t> massAssembler(geometry,basisDisplacement,bcInfo,gravity);
98 massAssembler.options().setReal("Density",density);
99
100 // creating time integrator
101 gsElTimeIntegrator<real_t> timeSolver(assembler,massAssembler);
102 timeSolver.options().setInt("Scheme",time_integration::implicit_nonlinear);
103 timeSolver.options().setInt("Verbosity",solver_verbosity::none);
104
105 //=============================================//
106 // Setting output & auxilary//
107 //=============================================//
108
109 // displacement field
110 gsMultiPatch<> displacement;
111 // stress field
112 gsPiecewiseFunction<> stresses;
113 // constructing an IGA field (geometry + solution)
114 gsField<> dispField(geometry,displacement);
115 gsField<> stressField(assembler.patches(),stresses,true);
116 std::map<std::string,const gsField<> *> fields;
117 fields["Displacement"] = &dispField;
118 fields["von Mises"] = &stressField;
119
120 std::ofstream logFile;
121 logFile.open("flappingBeam_CSM3.txt");
122 logFile << "# simTime dispAx dispAy compTime numIters\n";
123
124 gsProgressBar bar;
125 gsStopwatch iterClock, totalClock;
126
127 //=============================================//
128 // Initial conditions //
129 //=============================================//
130
131 // set initial conditions
132 timeSolver.setDisplacementVector(gsMatrix<>::Zero(assembler.numDofs(),1));
133 timeSolver.setVelocityVector(gsMatrix<>::Zero(assembler.numDofs(),1));
134
135 timeSolver.constructSolution(displacement);
136 assembler.constructCauchyStresses(displacement,stresses,stress_components::von_mises);
137 writeLog(logFile,displacement,0.,0.,0);
138 // plotting initial displacement
139 gsParaviewCollection collection("flappingBeam_CSM3");
140 if (numPlotPoints > 0)
141 gsWriteParaviewMultiPhysicsTimeStep(fields,"flappingBeam_CSM3",collection,0,numPlotPoints);
142
143 //=============================================//
144 // Solving //
145 //=============================================//
146
147 real_t simTime = 0.;
148 real_t numTimeStep = 0;
149 real_t compTime = 0.;
150
151 gsInfo << "Running the simulation...\n";
152 totalClock.restart();
153 while (simTime < timeSpan)
154 {
155 bar.display(simTime/timeSpan);
156 iterClock.restart();
157
158 timeSolver.makeTimeStep(timeStep);
159 timeSolver.constructSolution(displacement);
160
161 //assembler.constructCauchyStresses(displacement,stresses,stress_components::von_mises);
162
163 compTime += iterClock.stop();
164 simTime += timeStep;
165 numTimeStep++;
166
167 if (numPlotPoints > 0)
168 gsWriteParaviewMultiPhysicsTimeStep(fields,"flappingBeam_CSM3",collection,numTimeStep,numPlotPoints);
169 writeLog(logFile,displacement,simTime,compTime,timeSolver.numberIterations());
170 }
171
172 //=============================================//
173 // Final touches //
174 //=============================================//
175
176 gsInfo << "Simulation time: " + secToHMS(compTime) << " (total time: " + secToHMS(totalClock.stop()) + ")\n";
177
178 if (numPlotPoints > 0)
179 {
180 collection.save();
181 gsInfo << "Open \"flappingBeam_CSM3.pvd\" in Paraview for visualization.\n";
182 }
183
184 logFile.close();
185 gsInfo << "Log file created in \"flappingBeam_CSM3.txt\".\n";
186
187 return 0;
188 }
189